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ABSTRACT 

The multi-level (multi-grid) adaptive technique is a gen- 
eral strategy of solving continuous problems by cycling-, 
between coarser and finer levels of discretization. It pro- 
vides very fast general solvers, together with adaptive, 
nearly optical discretization schemes. In the process, 
bcwndary layers are automatically either resolved or skipped, 
depending on a control function vhieh expresses the computa- 
tional goal. The global error decreases exponentially as a 
function of the overall computational work, in a uniform rate 
independent of the magnitude (c) of the singular-perturbation 
terms. The key are high-order uniformly stable difference 
equations, and uniformly smoothing relaxation schemes. 

1. INTRODUCTION 

The Multi-Level Adaptive Technique (MLAT) is a general 
numerical strategy for solving continuous problems such as 
differential and integral equations and functional minimiza- 
tion problems. It will be discussed here mainly in terms of the 
numerical solution of partial differential boundary-value prob- 
lems, with special emphasis on singular-perturbation problems. 


The work reported here was perform’d under NASA Contract No. 
NASl-14101 while the author was in resi eiice at ICASE, NASA 
Langley Research Center, Hampton, VA ''.j665. 


The usual approach is first to discretize the boundary- 
value problem in some preassigned manner ^c.g., finitc-ele- 
nent or finite-difference equations on a fixed grid), and then 
to submit the resulting discrete system to some numerical 
solution process. In MLAT, however, discretization and solu- 
tion processes are intermixed with, and greatly benefit from, 
each other. A sequence of uniform grids (or "levels"), with 
geometrically decreasing mesh-sizes, participates in the pro- 
cess. The cooperative solution process on these grids in- 
volves relaxation sweeps over each of them, coarse-grid-to- 
fine-grid interpolations of corrections and f ine-to-coarse 
transfers of residuals. This process has two important basic 
benefits. On one hand it acts as a very fast general solver 
of the discrete system of equations (including the equations 
on the finest grid). On the other hand it provides, in a 
natural way, a flexible, adaptive discretization. For 
convenience, we discuss these aspects one by one. 

2.1. The Fast Solver 

In Section 2 of this paper we portrait the multi-level 
process as a fast solver, i.e., regarding the coarser grids 
as nothing but auxiliaries for solving the fines t-grid equa- 
tions. The description has appeared before (Brandt ''1972), 
(1977a), (1977b)), but here we add more detailed examples of 
the solution process (Sec. 2.2), and emphasize some new 
Important aspects. In particular, the "f ine-to-coarse cor- 
rection function" (or the "local relative truncation error") 

JJ 

T, is discussed, together with some of its usages, 
h 

One usage, developed in collaboration with N. Dinar, is 
the so called -^-extrapolation. It amounts to a trivial 
additional operation in the multi-grid program, and costs 
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nttgligible aaount of extra computing time. But» aa shown In 
Sec. 2.3, It improves the solution, sometimes by very much. 
With It, at a total computational cost of 4 to 7 work- 
units (a unit being the work equivalent of one Gauss-Seldel 
relaxation sweep over the finest grid), a solution u^ Is 
always obtained which Is better (i.e., a closer approximation 
to the true differential solution U) than the exact solu- 
tion U" of the difference equations. Furthermore, In case 
some extra smoothness Is present, u^ will be some orders- 
of-magnitude better than U^; or, alternatively, it may be 
obtained at a much smaller computational cost. 

As other usages of the f Ine-to-coarse correction function 

T. we list, in Section 2.4, some very efficient methods for 

Q 

aaJcing nonlinear continuation (e.g., for bifurcation 
problems) , methods for optimal-control problems, ill-posed 
boundary-value problems and parabolic time-dependent problems , 
as well as fast solution methods that can operate with a 
limited computer storage. Also mentioned in Section 2.4 are 
new msaerlcal experiments, made in collaboration with 
N. Dinar, for the steady-state incompressible Navier-Stokes 
equations, including the singular-perturbation case of large 
Reynolds numbers. These, and many other experiment:; 
briefly referred to, clearly indicate that the above-men- 
tioned multi-grid efficiency (solution in just few work- 
units) is obtained for general elliptic and non-elliptic 
systems on general domains. 

The fast-solver aspect of the multi-level techniques was 
studied by various other workers, starting, perhaps, with 
the "group relaxation" of Southwell (1935) . See references 
in Brandt (1977a), and more recent references in Nicolaides 
(1978) and Hackbusch (1978b) . Most of this work is very 
theoretical. That is, rigorous asymptotic bounds are 
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derlvtd for the sulti-grid efficiency. The price of 
rigorosity, of course, is that the results, are far from 
realistic: The proofs hold only for extremo.ly small mesh 

sizes, and, even for those, the work estimates are orders- 
of -magnitude too large. (Cf. Section 10 of Brandt (1977a).) 
The rigorous estimates are too crude, in fact, to yield any 
useful Information; e.g., they cannot resolve the difference 
between more efficient and less efficient multi-grid 
processes. For singular-perturbation problems, especially, 
the rigorous proofs hold only for mesh-sizes which are 
microscopic compared with the practical ones. For this 
reason, and since the quantity we try to estimate here is 
actually nothing but the computer -time, which of course we 
know anyway (at least aposteriori) , a different type of 
theoretical studies are preferred by the present author. 
Briefly, discarding rigorosity, it is observed that the 
important multi-grid processes are of local nature, since 
long-range convergence is obtained by coarse-grid processes, 
which cost very little. One can therefore analyze the 
crucial aspects of multi-grid processes by employing a local 
mode (Fourier) analysis, ignoring for instance far bound- 
aries and changes in the (possibly nonlinear) equations . 
Experiments with various types of equations (see Dinar 
(1978) and also Poling (1978)) shows that this analysis 
(which is much simpler than the rigorous theorems) precisely 
predicts the multi-grid efficiency. It is therefore very 
useful in selecting efficient algorithms (see, e.g.. 

Appendix A in Brandt (1977a)), in understanding the numeri- 
cal results, and in debugging multi.-grid programs. It led, 
in fact, to the efficient algorithms mentioned above, which 
solve, for example, Navier-Stokes equations in about 7 work- 


units. Its mechanized use if mentioned in Section 2.4 and* 
for singular-perturbation problems* in Section 6. 

About the difficulties in implementing multi-grid proce- 
dures* and vhat to do about then, some general coionents are 
made toward the end of Section 2.4! 

1.2. Kon-Vniform, Adeptible Structures 

Section 3 of this paper (following and adding to Sections 
2 and 3 in Brandt (1977b)) discusses the special capability 
of the multi-level structure to create non-uniform, flexible 
discretization patterns, especially such patterns as required 
by singular-perturbation problems, where thin layers should 
sometimes be resolved, either near or away from boundaries. 
This capability is obtained by observing that the various 
grids (levels) need not all extend over the same domain. 

Finer levels may be confined to increasingly smaller sub- 
domains, so as to provide higher resolution only where 
desired. Moreover, we may attach to each of these localized 
finer grids its own local system of coordinates, to fit curve 
boundaries or to approximate directions of interior inter- 
faces and thin layers. Unlike global coordinate transformation, 
these local coordinates do not complicate the difference 
equations throughout the domain (hence do not turn the one- 
dimensional trouble of boundary approxl^nations into a two- 
dimensional trouble of complicated equations). All these 
patches of local grids interact with each other through the 
multi-grid process, which, at the same time, provides fast solu- 
tions to their difference equations (an Important advantage 
over other methods of patching grids or using transformations) . 

This structure, in which non-uniform discretization is 
produced through the sequence of uniform grids, is highly 
flexible. Discretization parameters, such as (finest) 
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mcsh-sizes and approximation orders, can be locally adjusted 
in any desired pattern, expending negligible amounts of 
book-keeping work and storage. 

In particular, since in this structure only equidistance 
differencing is needed (much less expensive than differenc- 
ing on variable grids) , it becomes feasible to employ high- 
order difference approximations, even in singular-perturba- 
tion cases (see Section 5) . 

The discretization can thus be progressively refined and 
adapted. The actual adaptive solution process is governed 
by certain criteria, described in Section 3.6. Derived from 
optimization considerations, these are local criteria which 
automatically decide where and how to change the local dis- 
cretization parameters. Furthermore, these criteria are con- 
trolled by the user through a certain function G (the error- 
weighting function), which, in effect, expresses the purpose 
of the numerical calculations, i.e., the sense (or the error 
norm) in which approximations to the true solution are to be 
measured. The resulting discretization will be of high order 
wherever the evolving solution is suitably smooth. Singu- 
larities of all kinds will automatically be detected and 
treated, usually by introducing increasingly finer levels on 
increasingly smaller neighborhoods of the singularity. 

1.3. MLAT Solutions to Singular Perturbation Problems 

The discretization patterns produced by this general 
adaptive process for singularly-perturbation cases are 
studied in Section 4. It turns out that boundary layers 
are sometimes resolved by the adaptive process, and in other 
cases they are completely "skipped", depending on the choice 
of the control function G. The decision whether and how 
to resolve the layer is automatically taken by the adaptive 
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process itself. In any case, the convergence rate In the 
suitable sense (i.e., in the error norm corresponding to G) 
is always fast. 

Rates of convergence of adaptive processes are measured 
by the rate of decrease of the error £ ■■ ||u - u|| as a 
function of the computational work W, where u ■ u(W) 
is the evolving numerical solution, U is the true differ- 
ential solution, and 1| • 1| is the appropriate error norm. 
Since the grid is not uniform, nor does it have any fixed 
number of grid-points, the work-unit in this context must 
be different from the one mentioned above (which was defined 
in terms of the finest grid). It can be defined, for 
Instance, as the work of applying the lowest-order differ- 
ence equations at one grid-point. Thus, for example, in 
the conventional case, where regular grids are applied for 
solving b d-dimensional regular differential problem, 
employing O(h^) difference approximations, if a fast 
solver (e.g., a multi-grid algorithm) is used which solves 
the algebraic equations in 0(h arithmetic operations, 
then the rate of convergence can be expressed as 
E a The constant C depends not only on p but 

also on the solution U. 

We show that, using adf.ptive discretization, the same 
p-order convergence rate E ■ 0(W is obtained uni- 

formly in the size (t) of the singular perturbation; that 
is, C (p) does not depend on c. Moreover, if the order- 
of -approximation p is adaptiblc too, the rate of conver- 
gence is uniformly exponential. More precisely, for cases 
requiring boundary-layer resolution it is shown that 
£ cW^cxp(-cW°) , with a and c independent of e. We 
assume, of course, that the convergence rate for the reduced 
problem would not be slower. In cases the boundary-layer is 


skipped, the rate depends solely on the rate of the reduced 
problem. We also show that for this type of results it is 
not necessary to adapt p locally; p may be an adaptible 
constant. 

These convergence rates are uniform; they are obtained 
for any c, small, moderate or large. Nothing actually 
should be known in advance about the value of c. It is 
not even required at all to know that this is a singular- 
perturbation problem. Most other solution methods, by 
contrast, solve either the regular case (e “ 0(1)) or the 
asymptotic case (e very small), but not the whole range. 
(Quite often, intermediate values of r. are the most 
difficult to solve). Here, no special analyses are required, 
no need to separate the reduced problem from the singular- 
perturbation, and, In particulai', no need to compute the 
proper reduced boundary conditions. Mo matching procedures 
are employed. The method works similarly for interior 
singular layers, e.g., for ODE problems with turning points, 
even when (as in nonlinear problems) the location of the 
singularity is not known in advance. 

Although no apriori knowledge is needed about the size 
of the singular-perturbation, some rules should be kept in 
dealing with potentially singularly-perturbed problems. As 
is well known, difference schemes should be constructed 
which are uniformly stable. This aspect is discussed in 
Section 5, including the construction of high-order uni- 
formlg-stable difference operators. Similarly, in the multi- 
grid processing of potentially singular problems, relaxation 
schemes with uniform smoothing rates should be employed. 

Such schemes are described in Section 6. 

Remark. Sections 5 and 6 are abbreviated here. For 
their full text, see Brandt (1978b). It will include 
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further discussion of elliptlclty and uniform wcll-poscdness 
of flnlte>dif£ercnce systems, as well as remarks on the 
uniformly well-posed approximation of singular-perturbation 
problems with highly-oscillating solutions. 
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2. MULTI-GRID FAST SOLVERS 


To understand the basic numerical processes of MLAT, 
consider first the usual situation where a partial differ- 
ential problem 

LU(x) ■ F(x), for X • (x, ,...,Xj) in a domain (2.1a) 

A i Q 

n c R** , 

AU(x) ■ C(x), for X on the boundary of G (2»lb) 

is discretized in a preassigned manner on a given uniform 
grid . G^, with mesh-size h, yielding the finite-differ- 
ence equations 

lV'(x^) - F^'Cr^), (x*‘ c G**) . (2.2) 


Here U - (U, ) and its discrete approximation 

ij I - *1 

U^ are q-dimensional vectors of unknown functions, L and 
A are linear or nonlinear differencial operators and 

is, correspondingly, a linear or nonlinear expres- 
sion involving values of at x^ and at neighboring 

grid points. At various Instances of the solution process, 
we have on an approximation to U^, which we will 

generally denote by u^. 

In this section multi-level techniques for the fast solu- 
tion of (2.2), with coarser grids serv'ing as auxiliaries, 
will be described. In this context the term multi-grid, 
rather than multi-level, can be used. The difference 
between “grid" and "level" arises only in the more general 
situation (see Section 3 below). 
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2.2. Basic Halti-'Gr id Processes 


To obtaio a fast solution to equation (2.2) via the multi- 
grid method, VC add to a sequence of coarser uniform 

grids. Let C" be such a coarser grid; e.g., let the grid- 

1] V 

lines of G' be every other grid-line of C , so that its 
meshsite ij H - 2h. 

One way of inexpensively obtaining an approximate solu- 
tion to (2.2) is first to obtain an (approximate) solu- 

u 

tion u to the corresponding coarser problem 

l”uV) • r®(x^), (x®« G®) . (2.3) 

(which is much less expensive to solve since it contains far 
fewer unknowns) and then to interpolate u to the fine 
grid: 


h -H H 
u ■ I. u 

H 


(2.4) 


The symbol l" stands for the operation of interpolating 
H ^ h 

from G" to G . Polynomial interpolations of any order 
can be used. (The optimal order is discussed in Section A. 2 


of Brandt (1977a). 
of derivatives of 


Generally, if m is the highest order 


U 


j 


i 

in L and p is the order of 


approximation, then an interpolation of order at least 

m + p (i.e., polynomials of degree at least m, + p - 1) 

j H ^ 

should be used for to ensure full multi-grid efficiency. 

In some particular situations, even greater efficiency can 

be achieved by still higher interpolation; see footnotes 1 

and 5 below.) 

How good the approximation (2.4) is depends on the smooth- 
ness of the solution U^. In some cases U® is so smooth 
that, if the interpolation and the coarse grid operator 
L® are of order high enough to exploit that smoothness, 
then u® obtained by (2.4) saf**;fles 


-U- 


Ilu**- u|| < |ju^- U|| , (2.5) 

in some suitable norm. This means that solves (2.2) 

"to the level of the truncation error", which is all ve can 
meaningfully require in solving (2.2). In such cases, how- 
ever, the fine grid is not really needed: the coarser 

grid C already yields a solution with the required 
accuracy. If is at all needed, our first approximation 

(2.4) will require a considerable Isprovement. 

Can we compute a correction to u^ again by some inter- 

polation from the coarse grid C ? Namely, can we somehow 

h K h H 

approximate the error V ■ U ‘ - u by some v" computed 

HI) U 

on C ? Normally , the answer is no. If u in (2 4) is 

a good enough approximation to U^, then V will be a 

rapidly oscillating function chat cannot meaningfully be 

described on the coarse grid C . Therefore, before we can 

reuse coarse grids, the ctror should be smoothed out. 

An efficient smoothing is obtained by reI*x*tion sweeps. 

A standard example is the Causs-Seidel relaxation sweep. 

This is a process in which all points x^ of are 

scanned one by one in some prescribed order. At each point 

the old value u^(x^) is replaced by a new value, which is 

computed so that (2.2) is satisfied at that particular point 

(or nearly satisfied, in case (2.2) is nonlinear at chat 

point; one Ke^’ton step of changing u^(x) is enough). 

Having completed such a sweep, the system (2.2) is not yet 

solved, because its equations are coupled to each other; but 

the new approximation is hopefully "better" than the 

old one. 

In fact, a well known, and extensively used, method for 
solving (2.2) is by a long sequence of relaxation sweeps. 
When the system (2.2) is linear, convergence of u^ to U 
is obtained by a sequence of relaxation sweeps if and only 
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if the system is definite. But the rate of convergence is 
very slow. Typically, if m is the order of L and 
is the number of grid intervals in the direction, then 

the number of sweeps required for convergence is proportional 
to (min[N^,...,N^])“. 

A closer examination, e.g., by Fourier analysis, of the 
error V^, shows that the components slow to converge are 
those whose wavelength is large compared with the mesh-size 
of h. The high-frequency components, however, converge 
very fast; they practically disappear after a few sweeps. 

For example, in Gauss-Seidel relaxation for the 5-point 
Laplace operator 

L^**(x,y) = A^U^Cx.y) = {U^(x + h,y) + U^(x-h,y) 


+ U^(x,y + h) +U^(x,y-h) -4U*'(x,y)} , (2.6) 

the convergence factor of the Fourier component 
exp[i(6^x + 02 y)/h] (i.e., the factor by which the magni- 

tude of its amplitude in the .error expansion is reduced by 
one sweep) is 


v(ej,.e2) 


4 





(2.7) 


For the longest components S. = 0(N. ), and hence 
-2-2 J J 

p = 1 - + 1^2 ) . But for high-frequency components, 

say with max[e.[ i ? we have p ^ .5, so that in three 
relaxation sweeps these components are reduced by almost an 
order of magnitude. 

This means that relaxation sweeps, inefficient as they 
are in solving problems, are very efficient in smoothing out 
the error V^. This property, which is extensively used in 
multi-level algorithms, is very general. It holds for 
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Gauss-Scidel relaxation of any uniformly elliptic scalar 

(q “ 1) difference operator, whether linear or nonlinear. ^ 

For elliptic systems (q > 1), efficient smoothing is 
obtained by suitable variants of the Gauss-Seidel relaxa- 
tion. Even degenerate and singularly-perturbed elliptic 
operators are smoothed out with similar efficiency, provided 
more sophisticated variants are used, such as line relaxa- 
tions in suitable directions, or "distributed" Gauss-Scidel 
relaxations. Moreover, some of these variants are very 
efficient even for non-elliptic systems. (See Section 3 in 
Brandt (1976), Section 3 in Brandt (1977a), Lectures 5, 6 
and 7 in Brandt (1978a) and Section 6 below.) It is also 
important to note that fortunately the smoothing efficiency 
does not depend on some sensitive relaxation parameters.. 

Such parameters are sometimes needed (c.g., a relaxation 
factor is required in simultaneous-displacement relaxation 
schemes, which are used in conjunction with vector or 
parallel processing) ; but since smoothing is a local process, 
the optimal values of the parameters depend on the local 
operator only, and can easily be calculated by local Fourier 
analysis. Large deviations from the optimal values have 
only mild effect. 

; Thus, after a couple of relaxation sweeps, the error V 

! 

‘is smooth, and a good approximatipn to it can inexpensively 
’ ' H 

‘be computed on the coarser grid G . To see how this is 
done in the general nonlinear case, observe that on G^ the 
equation satisfied by is the "residual equation" 

I lN^(x^) «= r^(x^), (x^e G^) , - - (2.8a) 

^ where 

£hyh ^ + yhj _ j^h^h ^ <2. 8b) 

r^ S F^ - . (2.8c) 
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*h h 

In the Linear case L ^ L , and on first reading one nay 
like to keep this case in mind. (2.8) is of course fully 
equivalent to (2.2), but we are interested in this form 
because V^, not U^, is the smooth function which we like 
to approximate on the coarser grid G^. r is the 
"residual function", and, like V^, it is smoothed-out by 
relaxation. The approximation to (2.8) on the coarse grid is 


.H._H h . , H. h. , H. r:H h. H. 

L (I^u + V )(x ) - L (IjjU )(x ) = Ij^r (x ), 


(x” £ g”) , 


(2.9) 




where V is designed to be the coarse-grid approximation to 
h H "H 

V , and 1. and I. are interpolation operators (not nec- 

” ” h H h 

essarily the sane) from G to G . Since the points of G 

H 

are often a subset of the points of G , one can actually 

use direct '’injection", i.e., l||u^(x^) * u^(x^). In many 

cases, however, it is prcferrable to use "full weighting", 

H h H 

i.e., to use I u (x ) which is a weighted average of 
h,..h, h _h H 


values u (x ) at points x c G 
h, h. 


near 


in such a way 


that all values u (x ) equally contribute to the coarse- 
grid values. 

Observe that at this stage we could not approximate the 
equation L^(u^ + V^) = f^ on the coarse grid by the 
simpler approximation 

l“(iJu^ + = I^f^ , (2.9’) 

since the error of this approximation depends on the rapidly- 

oscillating part of u^, which may be large compared with 

the function we seek to approximate. In (2.8), by 

contrast, even if is small, the left-hand side is still 

approximately a linear operator in V^, and the left-hand 

2 ) 

side of (2,9) nicely approximates that linear operator. 

In fact, the coefficients of that quasi-linear operator on 


I 
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6 



the fine grid depend on u^, while on Che coarse grid they 

H h H 

have similar dependence on • Hence, if is a proper 

averaging operator, the coarse grid coefficients will auto- 
matically be a-'erages of the fine grid coefficients, so chat 
even if u is highly oscillatory, the coarse-grid equation 
is a proper "homogenization" of the fine-grid equation. (For 
discussions of homogenization see, e.g., Bnbuska (1973) and 
Spagnolo (1975).) 

Observe also that residuals are defined, and are 

transferred (with some averaging) to the coarser grid, not 
only with respect to the interior equations, but also with 
respect to the boundary conditions. In order that such 
transfers are done in the right scale, it is important that 

(i) the difference equations (2.2), and similarly (2.3), 
approximate (2.1) without change of scale (e.g., without 
multiplying through by h*. Equations (2.8)-(2.9) refer to 
the divided form of the difference approximations. Keeping 
this in mind, one can of course write the program with 
differently scaled equations, provided r^ is multiplied 

by a suitable constant when transferred to the coarser grid.) 

(ii) Difference-equations approximating different differen- 
tial equations should be clearly separated. For example, do 
not scramble together equations approximating (2.1a) with 
equations approximating (2.1b). Do not incorporate the 
boundary condition into the neighboring interior equation. 
Failure to observe these rules is a common error in multi- 
level programming. 

To avoid complicated linearization in solving (2.9), a 
new unknown function 

U? » + v” (2.10) 

n n 

is introduced (instead of V®) on G^, in terms of which 


K ■ 
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[ 


r 


(2.9) becootcs 

lXu") - (*” « o") . , (J.Ua) 

where 

fJ « L^(lju^) + . (2.11b) 

The advantage of this new form Is that it is the same equa- 
tion as (2.3), except for a different right-hand side. (The 
difference between the two right-hand sides is an important 
quantity which will be exploited below. See Section 2.3.) 
Moreover, (2.11) and (2.3) has the same form as (2.2). 

Bence, the same routines (e.g., the same relaxation routine) 
can be used in treating all of them. (See for example the 
simple sample program in Appendix B of Brandt (1977a) .) 

It is also worth noting that our new unknown (2.10) • 

represents, on the coarse grid, the sum of the basic approx- 

li li H 

imation u and its correction V . Thus, u. is the full 

H ^ 

current approximation, represented on G . The scheme of 

JJ 

using u, is therefore called the Full Approximation Scheme 
a 

(FAS) . To be distinguished from the Correction Scheme (CS) , 
which directly uses V^. (The Correction Scheme is messy in 
nonlinear problems, and cannot be applied on composite grids 
(see Section 3) . We therefore continue our discussion here 
only in terms of the more general scheme FAS.) At conver- 
gence, when i and « 0, we have 

Thus is a coarse-grid function which coincide with the 

h 

fine-grid solution - a fact which will also be very useful 

below (Sections 2.3 and 2.4). 

Once (2.11) is solved (or approximately solved), we want 

to use its solution to correct the basic approximation u^. 

^ H H h 

In doing so we should keep in mind that « U, - I. u , 

11 il 

and not u” itself, is the function which approximates a 
n 
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f 


I 


I 


smooch fine-grid function, and hence V (or its computed 
approximation) is vhat we should Interpolate to the fine 
grid and use it there as a correction to" u^. Thus, denot- 
ing our approximate solution to (2.11) by u^, the corrected 
approximation on the fine grid should be 


h 

“new 


“old ^H^“h 


h . 
^h“OLD^ 


( 2 . 12 ) 


Observe that this interpolation is not equivalent to 


-h H 

“new “ ^H“h ’ 


(2.13) 


Vi H 

since is not the identity operator. The important 

n n 

difference is that (2.12) preserves the high-frequency 

information contained in while (2.13) does not use 

“old' thus loses this information. Interpolation of 

the type (2.12) is called FAS interpolation. 

The order of the interpolation in (2.12) need not be 

as high as in the first coarse-to-fine interpolation (2.4). 

Order m. is enough (cf. the discussion following (2.4)). 

3 

For example, if the differential equation is of second-order, 

then linear interpolation is enough. 

We sunmarize the basic processes above: To solve the 

fine-grid eguations (2.2), an initial approximation (2.4) is 

H 

obtained from an approximate solution u to the coarse grid 
equation (2.3). Then the approximation is improved by a 
"multi-grid cycle”. This cycle includes a couple of relaxa- 
tion sweeps followed by the ”coarse-grid correction" (2.12), 
H 

in which uf~ is an approximate solution to the coarse-grid 
h 

correction eguations (2.11) . 

In most cases, at the end of the multi-grid cycle the 
approximation u^ will satisfy (2.5) and therefore require 
no further improvement. This is because the relaxation 
sweeps effectively liquidate the high-frequency 


I 
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components of the error, while the coarse-grid correction 
liquidates the lower components. In fact, we will see below 
(Section' 2.3) that, with some modification in the algorithm, 
at the end of une multi-grid cycle ||u - U|| may be much 
smaller than Hu^ - ul|. If, however, for any reason, a 
greater accuracy in solving (2.2) is desired, additional 
multi-grid cycles can be performed. Typically, each cycle 
which Includes three sweeps of a suitable relaxation scheme 
will reduce Hu^ - U^H by a factor of .2 to .08. 

We still have to specify how the coarse~grid equations , 

first (2.3) and later (2.11a), are actually solved. They 

are solved in the same way that (2.2) is solved, namely, by 

a combination of relaxation sweeps and coarse-qrid correct 

tlons, using a grid still coarser than G . More precisely, 

(2.3) is solved by a first approximation obtained from- a 

2H 

still-coarser grid (grid G , for example), and then a 

2H 

multi-grid improvement cycle (using G again) . For solv- 
ing equation (2.11a) the first approximation is l^u ; one 
multi-grid cycle (using G ) is enough for improving this 
approximation to the required level of accuracy. 

The full algorithm has several variations. One is flow 
charted here as Figure 1. This is essentially the same 
algorithm as described in Section 1.3 of Brandt (1977b). 
Sample runs of it can be produced by the MUGTAPE (1978a) 
program FASFMC. It contains three switching parameters: a, 

6 and n- Usually a * 2 where p is the approximation 
order. Optimal values for 6 and n are discussed in 
Sections A. 6 and A. 7 of Brandt (1977a). In practice the 
precise optimization is not important. One can take n * P 
and 6 = u , where v is an estimate for the smoothing 
factor (computed for example by SMORATE; cf. Section 6), 



'=‘ig. 1. FAS Full Multi-^rid (FAS FUG) Algorithm. 
(See Legend on next page.) 












Fig, 1. FAS Full Multi-Crid (FAS FMC) Algorithm. In this 
flowchart, as in the program itself, the different levels 

(grids) are not labelled by their mesh^size, as in the text, 
but by a positive integer k. k • 1 is the coarsest level, 
k “ M is the ultimately finest level, and k • i is the 
currently finest one (i.e., the finest so far used by the 
algorithm) . 

Thus, the original finite^difference equation on level k 
(before being changed to serve as a correction equation like 
(2.11)) is lV' - 7^. u* denotes the current approxima- 
tion, and f* th« current right-hand side, on level k. 

denotes interpolation of o^drr a (the order of the 
diffierential equation) from level k - 1 to the next fiver 

I 

level k. denotes higher-order (pi + p order) inter- 
polation. denotes transfer (by some averaging) -from 

level k + 1 to the next coarser level k. t, denotes an 

k 

approximate measure of the local truncation errors on level 

k (cf. Section 2.3). e^^ is a measure oi the residuals, 

taken during the relaxation sweep on level k. e^^ is the 

value of e. at the previous sweep, so that Cj. / e. > ti 

signals slow convergence, is a tolerance designed so 

that e^ < signals convergence of the current k-level 

problem. For k < Z the k-level problem is the correction 

problem to level k + 1, hence e, •» Se, On the current- 

k k+1 

ly finest level (k « Z) we need convergence to within the 
estimated size of the truncation error. Hence 
and before is known parameters 

H, 6 and a are discussed in the text. 


■B 




and r is the number of relaxation sweeps per multi-grid 
cycle. With good relaxation schemes p * .5 and r * 3. 

A slightly different algorithm is shown step-by-step in 
Table 1, in the next section. That algorithm can also be 
run through program FASFMG (using its FASFIX subroutine) . 

The difference between the two is that the algorithm in 
Fig. 1 is "accommodative'' t its flow depends on some internal 
checks, while that in Table 1 is "fixed", its entire flow is 
prescribed in advance, depending only on some input para- 
meters. 

2.2. Full Multi-Grid Run: An Exeunple 

Table 1 shows the steps and the results of a multi-grid 
solution process for the 5-point Poisson equation (cf..(2.6)) 

AJ^U^(x^ - F(x^, (x^e , (2.14) 

where is a 97 x 65 grid with mesh-size h * 1/32, 

covering the rectangle tO Ji 3, 0 _< X 2 £ 2} . Dirichlet 
boundary conditions are 'given on the boundary of the rectan- 
gle. In this particular run the boundary conditions and 
F “ AU were chosen so that the exact solution to the 
differential problem is U « sin(3Xj^ + 3 X 2 ). 

The flow of the algorithm can be seen from the first 
column of the table. It tells us the mesh-size H of the 
grid on which a Gauss-Seidel (GS) relaxation sweep is made. 
Thus, the process starts with 5 sweeps on G^^^, a 7x5 
grid with mesh-size H * 16h = 1/2, starting with the 
approximation u =0. Since the grid is very coarse, 
after 5 sweeps u^^^ solves (2.3) well below the trunca- 
tion level. This last u^^^ is then cubic- interpolated 

8h 

(i.e., interpolation of order 4) to the finer grid G to 

8h 

serve there as the first approximation u , as indicated 
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t;.ble I 


Output of Nultigrid Runs for a Poisson Problem, Produced 
by the MVGTAPP (2978) Program FASFMG 


Level H 


l|r"|L 


Relaxation 

Work 


c 

a: 

• 

•'ll. 

Usual 

T-extrap. 


l|u”- U|| 


2.11 

.750 

.270 

.104 

.0417 


Cubic interpolation 


.0039 

.0078 

.0117 

.0156 

.0195 


.25 


.25 


.249 


8h -8h 16h 

“ ‘ h6h“ 


8h 

3.73 

.0351 


8h 

1.34 

.0508 

.150 

16h 

.570 

.0547 


16h 

.269 

.0589 

.0611 

32h 

.203 

.0596 


32h 

.0262 

.0605 


32h 

.00164 

.0615 

.0586 

16h 

.145 

.0654 

.0765 

8h 

.492 

.0811 

.0798 


.150 


16h 

^8h 




Cubic interpolation 


4h ,4h 8h 
“ " ^8h“ 


.0600 


.0538 

.0600 

.0499 


.0572 


4h 

1.01 

.144 



4h 

.602 

.206 

.0645 

.0407 




8h , 8h,, 

8h 

.418 

.222 



8h 

.288 

.237 

.0398 

.0246 

16h 

.150 

.241 



16h 

.0473 

.245 

.0117 

.0517 

32h 

.0114 

.246 



32h 

.000095 

.247 











TABLE I - Continued 




7i 


32h 

.000006 

.248 

.00713 

.00458 

16h 

.0145 

.252 

.0123 

.00565 

8h 

.0988 

.266 

.0162 

.00731 

4h 

.201 

.330 

.O'iSO 

.00525 


Cubic Interpolation u 


2h 


^2h 

^h“ 


2h 

.267 

.580 



2h 

.205 

.830 

.0168 

.00467 




4h , 4h,. 




"2h * 

4h 

.177 

.893 



4h 

.158 

.955 

.0146 

.00387 

8h 

.109 

.971 



8h 

.0752 

.986 

.00893 

.00242 

16h 

.0378 

.990 



16h 

.0139 

'<94 

.00206 

.000870 

32h 

.00551 

.995 



32h 

.000715 

.996 



32h 

.000045 

.997 

.00223 

.000109 

16h 

.00435 

1.00 

.00254 

.000253 

8h. 

.0274 

1.02 

.00364 

.000226 

4h 

.0504 

1.08 

.00404 

.000675 

2h 

.0665 

1.33 

.00419 

.000472 


Cubic interpolation 

2n 


h 

.0680 

2.33 

• 


h 

.0559 

3.33 

.00412 

.00448 




2h , 2h„ 
\ * ‘’h 

2h 

.0540 

3.56 



2h 

.0513 

3.83 

.00394 

.000411 

4h 

.0455 

3.89 



4h 

.0393 

3.95 

.00340 

.000313 

8h 

.0278 

3.97 



8h 

.0192 

3.99 

.00206 

.000167 

ICh 

.00945 

3.99 



16h 

.00355 

3.99 

.000562 

.0000594 

32h 

.00138 

3.99 



32h 

.000123 

3.99 



32h 

.000008 

4.00 

.000540 

.0000231 


h, 


.0138 


.00344 
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TABLE I - Continued 


16h 

.00113 

4.00 

.000625 

.0000261 


8h 

.00684 

4.02 

.000879 

.0C00176 


Ah 

.0127 

4.08 

.000984 

.0000337 


2h 

.0163 

4.33 

.00102 

.0000438 


h 

.0182 

5.33 

.00103 

.0000518 

.00086 


8h 

in the table. Two CS sweeps are then nade on C . Then, 

the table shows, a switch is aade back to the coarser grid 

Such a "switch to the coarser grid" eeans that 

coarae-grid correction equations such as (2.11) arc set up 

on namely, their right-hand side is calculated and an 

initial approximation u^^^ ■ is introduced. Then, 

X6h 

as shown, 2 relaxation sweeps are made on chose C equa- 
tions starting with chat approximation. Then a switch 'is 

32& 

made to the still-coarser grid C (our coarsest grid 

here, which is 4 * 3 and has therefore only 2 interior 

points), where 3 sweeps are made. Then a switch back to the 

finer grid is made. Such a "switch to the finer grid" 

means that a FAS interpolation, like (2.12), is made from 
3 X 

C to C . These interpolations are all of order 2, 

that is, lineAT interpolations. One relaxation on 

8h 

then switch back to G and one sweep on it, and a multi- 

8h 8h 

grid cycle for G has been completed. At this point u 

solves (2.3) to the level of its truncation error, so we can 

already use it, with cubic interpolation, as the initial 

All 

approximation on G . Etc., etc., the first column in 
Table 1 shows the flow all the way to a final approximation 
on the finest grid (k •« h) . 

H H H 

At each relaxation jweep over G the residualr r (x ) 

are measured and their norm is accumulated. This norm is 

shown on the second column of Table 1. The precise 
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definition of the measured quantity |]r ]|q is not impor- 
tant at this point. In this specific run vc chose the fol- 

lowing definition: r are the dynamic residuals, i.e., 

H K 

r (x ) , as defined by (2.8c), is calculate:* using values of 

H H 

u as they stand imnedi.’>tcly before the point x is 

scanned by the relaxation sweep, that is, immediately before 
H H 

the value u (x ) is changed. This type of residual is 

least expensive to calculate, since it is (almost) calculateu 

H H 

anyhow in the course of computing the nev^ value of u (x ) . 
The norm ||*||q used in the table is the grid analog of the 
continuous nom 

,, „ / C(x) |r(x) |dx dx 

llr|L - — ; ^ (2.15) 

/ C(x)dXj^dx2 

where C is some function related to Che error functional 

vc are interested in (see Section 3.6). In this particular 

2 

run we had C(x) * C , where C is the distance of x 

X X 

from the nearest corner. 

The third column of Table 1 shows the accumulated work, 

measured in work units. A work unit is the work of one 

relaxation sweep on the finest grid C^. Hence, a relaxation 
NTi ^2 

sweep on C is counted as N work units, since it con- 

^2 

tains about N the number of grid-points contained in G . 
This work count neglects any other work except that of relaxa- 
tion sweeps, because (i) Relaxation sweeps account for at 
least 70S of the total actual work, (ii) The other work, 
mostly that of interpolations, is not directly expressible 
in terms of work units. We could measuic running time, but 
this depends too much on the particular hardware, software 
and program being used, (iii) The theoretical prediction of 
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convergence races Is aadc*. too, in Ceras of these work units, 
so it is convenient to use then for comparing experiments 
with theory. 

Observe how little is the work accumulated on the coarser 

grids. In sum, this solution algorithm requires only 3.33 

work units. A more precise count of all operations (includ- 

H 3) 

ing interpolations) shows that, if jjr j] is not computed 
this algorithm requires less than 42n arithmetic operations, 
where n is the total number of points in the finest grid. 

Incidentally, in this particular problem (2.1A), most of 
these arithmetic operations are additions. In fact, one can 
arrange it so that the only multiplications (and divisions) 
involved arc by factors 4, 2, 1/2 or 1/4, which in binary 
floating-point arithmetic can be performed as additions. 

(For this purpose, cubic interpolation should be performed 
through the difference operator itself, as for example in 
Hyman (1977 j ,'j Thus, no multiplications or divisions are 
required. 

In experiments at ICASF made by Craig T. Poling, the time 
measured for this algorithm on a CDC 6600 was .083 seconds 
for a 33 33 grid, and .303 for a 65 65 grid. A similar 

algorithm (using another kind of relaxation) on CDC STAR-100 
computer required .011 seconds for a 65 65 grid and .0347 

seconds for a 129 129 grid, i.e., about 2 microseconds per 

grid point. 

The first three colianns in Table 1 exhibit the standard 

output of a multi-grid run. The last three columns can be 

produced only in experimental runs made for cases in vdiich 

the €ixact differential solution U is known. The purpose 

of these columns is to show the quality of the computed 
H H 

approximations u . Here, u denotes the current approxi- 
H 

mation on G . Thus, In the coarse of a correction cycle 
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I 



U 

for level h, < i', u denotes the full approximation on 
H 

level G , similar to (2.10). The exact solution of the 

H H 

difference equations (2,3) on grid G is denoted by U . 

The fourth column of Table 1 shows, at various stages, 

the maximum difference between the exact differential solu- 

H H 

tioii U and the c imr.Mted u , where the values of u at 
each stage are those obtained at the end of the G*^ relaxa- 
tion sweep corresponding to that row of the table. The 
fifth column of the table will be explained in Section 2.3. 
The sixth column gives, for each level H, the maximum 
discretization error |U - u|. It is shown on the row 
corresponding to the end of the multi-grid correction cycle 
for level H. The main observation is, of course, that at 
this stage j|u^- u|[ is not considerably bigger that 
||U - Ujj, so the algorithm indeed solves the discrete, 
problem to the level of the discretization error. 

This performance is predictable, and will be the same for 
any other data. The local mode analysis shows that the multi- 
grid cycle used here reduces ||u^ - U^|j by a factor .08. 

A factor .25 would in fact suffice, since all we need in this 

2li 

cycle is to reduce the error from approximately [lu - U[| 
to approximately ||u^ “ Ull* 

Moreover, observe the row one before the last in Table 1. 

2h 

It shows that already at this stage u solves the problem 
to the level of the h discretization error, i.e., 

||u^^- ull is not much bigger than u||. Hence, the 

level of the h discretization error is actually obtained 
in only 4.3 work-units. The total number of arithmetic 
operations (additions and subtractions) required is only 
about 35n. The full 42n operations are required only if 

we need the solution, to the same accuracy, at all points 

h 2h 

of G", not only at points of G 




- 28 - 



1 


Another interesting comment which follows concerns the 

storage requirement for this algorithm. If the algorithm 

indeed terminates at 4.33 work-units, no FAS interpolation 

to is madi;. Hence the values of need not be stored. 

All the operations made on are the initial cubic inter- 

polation from G , followed by two relaxation sweeps and 

2h 

the calculation of F. . All these operations can be made 

by one pass over G , requiring in storage only 5 columns 

of the grid at a time. (The first sweep of relaxation is 

made on the fourth of these columns, then the second sweep 

can already be made for the third column and the residual 

transfer can then be made for the second colinnn.) Thus, the 

algorithm requires no storage (not even external storage) 

for the finest grid. The storage required for the coarser 

grids is only y + ••• < that of the finest grid. A 
4 xD 3 

modified multi-grid algorithm can work with even smaller 
storage (see Section 2.4). 


2.3. Relative Local Truncation Errors and t Extrapolation 

To realize further aspects of the multi-grid method, we 
now slightly shift our point of view. Going back to the 
coarse-grid correction equations (2.11), we rewrite them, 
using (2.8c), in the form 


lV = + r? 


(2.16) 


where 
H 


fH ^ ^Hj.h 

h 


and 


H h. yh. h h. 

Th “ L (Ij^u ) - I^(L u ) . 

At convergence, where u^ = U^, 


(2.10), U 






'h ^ 
local truncation error of 


H 


(2.17) 


we have V =0 and, by 


Hence, r actually represents the 


,H 


relative to G 


i.e., the 


error which arises when the fine-grid solution is 
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substituted in the coarse-grid equation (2.3). Compare this 
to the usual local truncation error, i.e., to the error T 
which arises when the true differential solution U is 
substituted into the G equations (2.3); namely 

- L% - LU . (2.18) 

We can now reverse our viewpoint. Instead of regarding 
the coarse-grid as a devise for calculating the correction 
(2.12) to the fine grid solution, we can view^the fine grid 
as a devise for calculating the correction to the 

coarse-grid equations, a correction which will make the 
solution of these coarse-grid equations to coincide (up to 
interpolation) with the fine-grid solution, x^ is there- 
fore called the fine-to-coarse correction function. It is a 
of defect correction (cf. Sec. 3.4). This point of 
view open many more algorithmic possibilities, such as the 
multi-grid method for non-uniform discretization (Section 
3.3 below), continuation techniques, methods for ill-posed, 
optimal-control, and time-dependent problems, and small- 

storage procedures (Section 2.4). 

The quantity x^ itself is very useful. First, it is an 
approximation to the local truncation error x . More pre- 
cisely, we see from (2.17)-(2.18) that 

(2.19) 


H H ^h 

X, X - X 

n 


This relation will be used in the adaptive processes (Sec- 
tion 3.6). Here we show how to use it to improve the multi- 
grid solution. 

The local truncation error can be expanded in Taylor 
series, yielding always a relation of the form 


x^(x) = C(x)h^ + O(h^), (P > P) . 


( 2 . 20 ) 


where C(x) depends on the (unknown) solution, but does not 
depend on b. Applying (2.20) also to x^. and using 
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hence 


(2.19), we find 

H _ 1 H . « ...p. . h 

^ “ T. + 0(H*^), = « t 

1 - ”, 

usually n = "t) . (2.21) 

H ^ 

Thus, if ve replace t. in the coarse-grid correction equa- 
H ^ p 

tions (2.16) by t"/( 1 - n*^) . we effectively introduce the 
n • 

true truncation error (up to order H^) instead of the 

relative one, and the solution U? should become a much 

h 

better approximation to the true differential solution U. 

We call this replacement a local truncation extrapolation^ 

or, briefly, -extrapolation. 

Note that t, is defined with respect to both the inte- 
n 

rior difference equations and the boundary conditions, hence 
T-cxtrapolation can be applied for both. 

The T-extrapolation costs very little. Only one operation 
(multiplication by (1 - q^) ^) is added, and only at 
coarser - grid points, so it amounts to less than -j opera- 
tion per fine-grid -point. The stages in the algorithm where 
this extrapolation is made are shown in the fifth column of 
Table 1. Also shown in that column are the results of this 
operation in terms of the error []u^ “ uj] at various 
stages. We see that with exactly the same flow , the 
algorithm produces now much smaller errors; in fact u^ now, 
after 4.3 work-units, is a much better approximation (to U; 
than the exact finite-difference solution U^. 

W’ith more work and some changes in the algorithm (quintic 
Instead of cubic interpolations, and two instead of one 
multi-grid cycles for each level, T-extrapolation being made 
on the second of the two), ||u^ ~ Uj] could be made much 
smaller yet^^ . 

The impressive improvement depends of course on the 
smoothness of the solution. Similar improvements could be 
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obtained, in some such eases, by using Richardson extrapola~ 

tion (extrapolating from the solutions and u^) . Indeed, 

in case the solution is not smooth, e.g., if it oscillates 

vildly on the scale of G^, the x-extrapolation does not 

H 

considerably improve u . But exactly in this case the one 

multi-grid cycle is enough to reduce Hu^ - U^H well below 

||U - U|I, since, exactly in this case, t is not consider- 

H 

ably smaller than x . So the point of the x-extrapolation 
is that it can always be used, for negligible extra cost in 
either programming or computer time, and it produces u^ 
that, at 4. 3-5. 3 work-units, is guaranteed to be no worse 
than U^, the full solution of the difference equations, 

W’ith the nice additional feature that any available smooth- 
ness is automatically exploited to improve u^ even further. 

It should also be pointed out thr.t the x-extrapolation 
can improve the solution in many cases where the Richards-^n 
extrapolation cannot. The x-extrapolation depends on Taylor 
expansion for the local truncation error, while the 
Richardson extrapolation requires such an expansion for the 
global discretization error - U. In many cases the 
later expansion does not exist; for example, no such expan- 
sion is possible when the local truncation error is not 
uniform. Another nice feature of the x-extrapolation proce- 
dure is that it produces the improved solution at all points 
of the fine-grid, while Richardson extrapolation gives it 
only at points common to the fine and the coarse grid. 

A remark on both types of extrapolations: When the extra- 

polated solution is considerably better than the fine-grid 
solution, its accuracy is actually comparable to the accuracy 
of a higher-order solution on the coarse grid. (Because the 
coarse-grid higher-order error term is not removed by the 
extrapolation.) That higher-order coarse-grid solution is 
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In principle less expensive in computer resources than the 
extrapolated solution, since it docs not use the fine grid 
at all. A fully adaptive procedure (cf. Section 3.6) would 
probably prefer in such a situation to use the higher-order 
scheme on the coarser grid. 

2.4. Cenezal Properties of Multi-Grid Solutions (Advertise- 
ment) 

We have discussed at length the multi-grid solution of one 
particular problem. The algorithm, however, does not depend 
on any of the particular features of that problem. 

The shape of the domain is immaterial. Only low-frequency 
error components are affected by it, and such components are 
always liquidated on the coarsest grids for negligible 
amount of work. The only effect of complicated boundaries 
is to complicate the difference equations at adjacent points 
and thus to make each relaxation sweep somewhat more expen- 
sive. In terms of work-units as defined above, the effi- 
ciency remains the same. Many experiments (reported in 
Brandt (1972), with many more examples in Shiftan (1972)) 
confirm this. 

Variations in the coefficients, or nonlinearities , in the 
differential equations usually also affect only low-frequency 
components, and are therefore still treated at the same 
multi-grid efficiency. When these variations are wild, i.e., 
when the coefficients change significantly over the scale of 
the grid, attention should be given to the proper choice of 
residual-weighting (l!| in (2.9) or (2.11)), since the 
residual function r** is considerably less smooth than the 
error V^. (See discussions in Section A. A of Brandt (1977a). 
More precise rules of residual-weighting a:> given in 
Lecture 17 of Brandt (1978a). The weights near boundaries 


LI 
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depend on the type of boundary conditions.) It Is also 
important for such cases that the coefficients of the coarse- 
grid difference equations represent local averages of the fine- 
grid coeff icie.'its. This, in fact, is automatically obtained 
if (for variable-coefficients linear problems) the difference 
equations on each grid are derived by suitably averaging the 
differential equations, or if (for nonlinear problems, where 
the coefficients depend on the solution) the solution- 
weighting (l|| in (2.9) or (2.11)) is a full weighting 
operator. 

Numerical experiments (Brandt (1978a) p. 17-7, and more 
details in Ophir (1978)) were conducted with difference 
equations of the form 


= 

h 

Cl. 6 

0,6 ■ 

*o6 


+ 

«h 

0,6* 

-1 - 2U1 

^ ^06 


h 


, - 2U^ - + I 

I 0.6 


c+l . B 


o.B ’ 


( 2 . 22 ) 


where the coefficients a^ and c^ vary wildly; e.g., a^ 
and c^ are random; or a^. “ * 1 st even points 

Op Op 

(o + 6 even) but a * .01 and c = 1 otherw’ise; etc. When, 
and only when, the algorithm (with proper line relaxation) 
used full weighting in the transfer to the coarse grid of 
loth the residuals and the coefficients, the solution was 
almost as efficient as in corresponding constant-coefficient 
cases. 

Many experiments were also made for various nonlinear 
problems. In 1975-76 Multi-grid programs were developed for 
the steady-state small-disturbance transonic flow problems 
(see South and Brandt (1976), and Section 6.5 in Brandt 
(1977a)), in which the differential equations are mixed 
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elliptic-hyporijoJic, and the solutions contain shock dis- 
continuities. Although these programs are somewhat obsolete 
(with the present stage of multi-grid experience and know- 
ledge they could be improved in various ways) , they do 
clearly show that the typical multi-grid efficiency can be 
obtained for this type of problems. , 

Recently, multi-grid codes for steady-state incompressible 
Navier-Stokes problems have been developed (see Lecture 7 in 
Brandt (1978a), Brandt (1978c) and Dinar (1978)). The 
algorithm is similar to the one shown in Section 2.2 above, 
except for the more elaborated "distributed" relaxation 
scheme which is required here (a) because it is an elliptic 
system, not a scalar equation, and (b) because, for large 
Reynolds numbers R the system is singularly perturbed. 
Cavity and pipe problems were solved , with R ranging from 
0 to large values (but below the values causing instabili- 
ties). For any such R, the process required 6 to 7 work- 
units, and produced solutions closer to the true differential 
solutions (in cases those were known) than the exact solu- 
tions of the finite-difference equations. 

Successful multi-grid applications are also reported for 
the minimal-surface equations (D. J. Jones, Lecture 15 in 
Brandt (1978a)), for the pressure iteration in Eulerian and 
Lagrangian bydrobynamics (Brandt, Dendi and Ruppel (1978)), 
and for some simple problems in finite-element formulations 
(Nicolaides (1978) and Poling (1978)). Not listed here are 
some other multi-grid applications, which seem not to have 
realized the true multi-grid efficiency^\ 

In some nonlinear problems a continuation (or embedding) 
process should be made, either because there are several 
solutions to the differential equations or because an 
initial approximation to the solution cannot otherwise be 
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¥ ^ 

I obtained. In such cases a certain problem parameter, y 

I say, is introduced, so that instead of a single isolated 

I 

\ problem ve confiider a continuum of problems, one problem 

P(y) for each value of y in an interval Yq i Y i Y*» 
f where P(Yq) is easily solvable (e.g., it is linear) and 

P(y^) is the target (the given) problem. The continuation 

1 

I method is to advance y from Yq to Y* in steps 6 y* At 

I each step vc use the final solution of the previous step 

P(y ■ ^y) (or extrapolation from several previous steps) as 

an initial approximation in an iterative process for solving 

P(y). The main purpose is to ensure that the approximations 

ve use are all ”close enough” to the respective solution, so 

I that some desirable properties are maintained and -convergence 

is ensured. The process should of course be made carefully 

in case of bifurcation, when several solutions branch off at 

a certain value of y. (See, e.g., Keller (1977).) 

With the multi-grid solution method the continuation 

process may cost very little. First, because it can be made 

on coarser grids. Sometimes, however, too coarse grids do 

not retain enough of the solution features, and the continua- 

iL tion may not accomplish its purpose. This means that oscil- 

S lations on the scale of a certain mesh-size h cannot be 

ignored. These oscillations, however, do not usually change 

I much at each 6 y step. The trick then is to use form (2.16) 

of the coarse-grid equations. We can make several continua- 

H 

tion steps on the coarse grid G only, retaining the values 
H 

of r. fixed. By doing this we effectively freeze the high- 
i h 

[ frequency part of the solution, and retain its influence on 

I the coarse-grid equations. Only once in several (sometimes 

; many) steps ve need to calculate on the finer grid too, in 

H 

order to update the fine-to-coarse correction r. . This of 

h ” 

course may be carried further: the G equations themselves 
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may Include a tf’,- correction term, and once in several 
visits to G" we may go to calculate on C ' in order to 
update that term; etc. 

A similar technique can be used in optimal control 
problems. Here, some parameters controlling the partial- 
differential problem should iteratively be adjusted so as to 
achieve some optimal condition (minimal "cost") . In such 
problems ve can again make most of the iterations on coarse 
grids, with frozen values of the fine-to-coarse corrections 

making only infrequent visits to finer grids for up- 
h H H 

dating T^. (Together with we should also freeze the 

fine-to-coarse correction of the cost functional.) If, for 
example, the control Itself is a grid function (e.g., a term 
in the right-hand side of the equations) , a change in the 
control value at a point will introduce only smooth changes 
in the solution in regions away from that point. It is 
therefore enough to keep refined only a small portion of the 
domain at a time, while at other regions the coarse level is 
used with frozen fine-grid corrections. This technique 
should combine with the usual multi-grid approach of obtain- 
ing the first approximation on each finer grid by interpolat- 
ing from a solution to a similar problem on the next coarser 
grid. 

Such techniques can also serve some ill-posed problems , 
where the solution should fit some data which are not the 
normal, well-posed boundary conditions. In such situations 
only smooth components of the solution can be meaningfully 
fitted to the data (or to smooth averages of data) . The 
data-fitting can therefore be made on a coarse grid, G 
say, but the coarse-grid equations should have the fine-to- 
coarse correction (2.17), so that they represent a reason- 
able approximation to the differential equations. 
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SmalJ-storago multi-grid tlgorithms are also based on the 
fom (2.16) of the coarse-grid equations. Observe that, 

U 

once T is known, the fine grid is no longer needed. But 

H " h 

T, depends mainly on high-frequency components of u , 
n 

which can be computed by having in storage only a small seg- 
ment of G^. An algoritlkm based on this idea ("segmental 
refinement". Section 7.5 in Brandt (1977a)) requires in 
principle only some IS^log n storage locations, where d 
is the problem dimension and n is the number of points in 
the finest grid. No external memory is assumed. 

rime-dependent problems, especially of parabolic type, 
usually require the use of implicit difference equations. 

The system of equations to be solved at each time step is 
similar to the steady-state equations, and can be solved 
usually by one multi-grid cycle, starting from the previous- 
time solution as the first approximation. Moreover, in some 

important cases (e.g., the heat equation) after a short 
2 

time tjj ■ 0(h ) the high-frequency components 
of the solution practically reach their steady 
state, and further changes occur only in the smooth compo- 

u 

nents. Hence, for many time-steps, the values of t, hardly 

H ” 

change. Freezing for several time-steps we can then 

solve the coarse-grid equation (2.16) without using the fine 

grid at all. Only infrequently should we make a time-step 

with fine grid calculations, to update the high-frequencies 

H 

of the solution and the values of In this way we retain 

fine-grid accuracy but each implicit time step of this kind 
costs on the average much less than an explicit time step on 
the finest grid. This kind of techniques for parabolic tine- 
dependent problems, are discussed in Lecture 9 of Brandt 
(1978a), and some preliminary experiments are reported by 
Dinar (1978). 
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Parallel or vector processing can be fully exploited by 
the multi-grid algorithm. Tlic main processes, namely, inter- 
polations and relaxation sweeps, are completely paralleliz- 
ablc, although it requires the use of a slightly different 
type of relaxation, with smoothing rates (per sweep) some- 
what slower than the (sequential) Causs-Seidel relaxation. 
(See Section 3.3 in Brandt (1977a).) 

Implementation difficulties. A fair amount of knowledge 
is required in implementing multi-grid algorithms, including 
some general knowledge common to all multi-grid programs, 
plus particular expertise related to the specific type of 
problem at hand. 

Generally, one has to be familiar with the basic rules of 
interpolation and residual-weightl.'S, with the normal flow 
of multi-grid runs^\ and, last bur aot least, with the data 
structure used in multi-grid programs. Vithout a suitable 
data structure the program will become complicated, wich 
many unnecessary repetitions. It is advisable to follow the 
programming techniques exhibited in the simple Sample Program 
(Appendix B of Brandt (1977a)) and in the progr.’uas of 
MUCTAPE (1978). With this technique, each operation (such 
as relaxation, coarse-to-f ine interpolation, f ine-to-coarse 
residual weighting) is written once for all grids. Moreover, 
most operations can be written once for all programs; that 
is, the same code can be used by all the programs which use 
the same data structure. This includes the coarse-to-f ine 
and f Ine-to-coarse transfer operations (e.g., interpolation), 
the operation of introducing the values of a given function 
into a given grid, operations of generating and manipulating 
grids (e.g., augmenting a grid, coarsening a grid, generat- 
ing the Interior part of a grid, or transposing a grid from 
row structure to column structure), displays of grids or 
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grid-functions, etc. Three different data structures are 
used in existing programs: One for rectangular grids, cne 

for single-string grids (where the domain can be defined as 
{(x,y)Ixj^ ± Xji fj^(x) ±y ± fjCx)}), a .d one for general 
grids. The latter is called Che QUAD structure, and is 
described in Brandt (1977b) and in Lectures 12 and 13 of 
Brandt (1975a). With these techniques, the progranning of a 
miltJ~gri<3 solution for c new problem is essentially reduced 
to the programming of the relaxation routine, (The residual- 
weighting routine should also be programmed anew for each 
problem, but its part that depends on the problem is a 
simple modification of a similar part in the relaxation 
routine.) 

Particular expertise is required in designing the relaxa- 
tion sweeps. For a uniformly elliptic scalar equations the 
simplest Causs-Seidel is the best scheme (on sequential 
machines), but suitable modifications are required for 
degenerate-elliptic, non-elliptic, indefinite or singular- 
perturbation problems, as veil as for cases of parallel or 
vector processing, and for problems involving more than one 
unknown function. Generally speaking, the particular know- 
ledge required for designing relaxation is similar in each 
case to the specific knowledge required in discretizing the 
differential equations. Similar - but not identical. As ii. 
learning discretization methods, one should learn relaxation 
methods gradually, starting from simplest models, gaining 
some bfisic insights, and only then proceeding to compi x 
real-world problems. 

The design of relaxation is much facilitated by a standard 
gauge we have for apriori measuring the relaxation efficiency. 
The only role of relaxation in multi-grid programs is to 
smooth the errors. The efficiency of the entire algorithm 
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depends on (and can be approximately predicted by) the 
"smoothing factor" of the relaxaficn sweeps. Since smooth- 
ing is a local process, the smoothing factor can be calcu- 
lated by a loc?l mode analysis (cf. Section 6). For simple 
equations this can be done by hand (see examples in Sections 
?.l and 6.2 of Brandt (1977a), and Section 3 of Brandt 
(1976).). For general equations it is done by using the 
computer routine SNORATE, available on ^nJGTAPE (1978) . The 
user of this routine supplies a description of the relaxa- 
tion scheme and other parameters (in a format explained by 
comment .a ds in the routine) . The output contains the 
smoothing fact' r and other useful information, including an 
estimate of the multi-grid convergence factor per work-unit. 
The routine can therefore be used to optimize relaxation, 
i.e., to select the best relaxation type and parameters from 
a given set of possibilities. 

Some general orientation concerning the relaxation of 
singular-perturbation problems is given in Section 6 below. 

A major advantage of the multi-grid solution process, in 
particular for singular-pertu'-bation and other irregular 
problems, is its full compatibility with adaptive processes. 
The reason is that the multi-grid process in itself is adap- 
tive: in adaptive processes mesh sizes are adapted to the 

computed solution; the multi-grid process goes one 
step further and employ mesh sizes adapted to the error in 
the computed solution. Let us now turn to these adaptive 
aspects of the multi-level techniques. 
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3. NON-UNIFOWI AND ADAPTIVE DISCRETIZATIONS 
3.1. Non-Uniformity Organized by Uniform Levels 



In principl,-^, the multi-grid solution process described 
in Section 2 could work equally well when the finest grid 

is non-uniform and even non-rectangular , with grid points 
at arbitrary locations. A relaxation with good smoothing 
rates on a general grid is obtained by employing all line, 
and marching, directions. The main difficulty with general 
grids, however, is practical: Merely to formulate and use 

the difference equations, let alone solve them, is compli- 
cated. It requires lengthy calculations and large memories 
for storing geometrical information, such as the location 
of each grid point, its neighbors, the coefficients of its 
difference equation, etc. The multi-grid processing of -such 
arbitrary grids generates additional practical difficulties 
since it requires the introduction of coarser grids and the 
grouping of grid points in grid lines (for line relaxation). 

These arbitrary general grids, hov?ever, are not really 
needed. We will show below a method of organization vjhich 
is less general but in which any desired refinement pattern 
can still be obtained, and easily changed, with negligible 
book-keeping and with difference equations always defined on 
equi-distant points. This flexible organization will 
naturally lend itself to multi-grid processing and to local 
transformations, and will lay the groundw'ork for efficient 
adaptation. 

It is proposed to organize non-uniform grids as "composite 
grids". A composite grid is a union of uniform subgrids 


./ih _2h _h _h/2 


(3.1) 


where the superscripts denote the mesh-size. The grids are 
usually positioned so that every other grid-line of is 
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a grid line of 


Unlike the description in Section 2, 


however, the subgrids are not necessarily extended over the 

entire domain fi: The domain of may be only a proper part 

,2h 


of the domain of G 


so that different degrees of refine- 


ment can be created at different subdomains. See Figure 2. 
Each G^ is extended, as a rule, over those subdomains 
where the desired mesh-size is roughly 1.5h or less. G^ 


may be thus disconnected, but its domain is always a sub- 
2h 

' The effective mesh-size at each neighbor- 


domain of G* 

hood will be that of the finest grid covering that neighbor- 
hood. Clearly, any desired mesh-size h can be approximated 
by some effective mesh-size h', where 0.75h ^ h* 1.5h. 
Mesh-sizes never require better approximation. 

The composite grid is very flexible, since local grid 
refinement (or coarsening) is done in terms of extending (or 
contracting) uniform subgrids, which is relatively easy and 
inexpensive to implement. A scheme for constructing, extend- 
ing and contracting uniform grids, together with various 
service routines for such grids (efficient sweeping aids, 
interpolations, displays, etc.) is described in Brandt 
(1977b) and is partly available on MUGTAPE (1978a) . One of 
its advantages is the efficient storage. The amount of 
logical information (pointers) for describing a uniform grid 
is proportional to the number of strings of points, and is 
therefore usually small compared with the number of points 
on the grid. Similarly, the amount of logical operations for 
sweeping over a grid is only proportional to the number of 
strings. Changing a grid is inexpensive, too. 

Moreover, this composite structure will at the same time 
provide a very efficient solution process to its difference 


equations, by using its levels (3.1) also as the multi-grid 




sequence (as in Section 2) . Each G will automatically 
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Fig. 2. A piece of non-uniform grid {A) and the uniform 
subgrids it is made of (B, C, D, E) . 
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play the role of the correcting coarse-grid whenever the 
h/2 

finer subgrid G is present (sec Section 3.3). 

In adaptive procedures, the sequence of subgrids will be 
kept open-ended, so that we can add increasingly finer or 
coarser levels, as needed. (Increasingly coarser levels may 
be needed if the problem’s domain G is unbounded and the 
bounded computational domain is chosen adaptively) . 

The coarsest subgrid should of course be kept coarse 
enough to have its system of difference equations relatively 
inexpensive to solve. Hence, there will usually be several 
coarse subgrids extending over the entire domain Si. That 
is, they will not serve to produce different levels of 
refinement, but they are kept in the system for serving in 
its multi-grid processing. 

There seems to be certain waste in the proposed system, 

because one function value may be stored several times when 

Ic 

its geometrical point belongs to several subgrids G . This 
is not really the case; The extra values are exactly those 
needed for the multi-grid processing. In the process, the 
different subgrids may have different values at the same 
geometrical point. Moreover, it is only a small fraction 
(2 *^) of the points that are actually being repeated. 

3.2. Vnisotxopic Refinements 

For singular-perturbation and other problems, it is some- 
times desired to have a grid which resolves a certain thin 
layer, such as a boundary layer. Very fine mesh-sizes are 
then needed in one direction, namely, across the layer, to 
resolve its thin width. Even when the required mesh-size 
is extremely small, not many grid points are needed, since 
the layer is correspondingly extremely thin. (See Section 4.) 
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Provided, of course, that the fine mesh-size is used only in 
that one direction. We need therefore a structure for mesh- 
sizes which get finer in one direction only. 

In case the thin layer is along coordinate lines (x «= 
constant) , we resolve it again by using a sequence of 
uniform grids, except that their meshes are no longer square; 
h^, the mesh-size in the x^ direction, may be very differ- 
ent from h^. See Figure 3. We still require all mesh-sizes 
to be binary multiples of some basic size h^, that is, on 
the k-th subgrid the mesh-size in the x. direction is 
k 

h^ ■ 2 Qq! integer, i «= l,...,d) . (3.2) 

For the multi-grid processing we require that for each such 
subgrid k, except for the very coarsest (k ■= 1) , there 
exists in the scheme a coarser subgrid, number £ ■= £(-k) 
say, such that for each 1 < i < d either n,„ *= n., or 
n 


i£ 


n., + 1, and 1 n.„ > I n., . Grid £ will be the 
xk x£ ik 


grid from w’hich corrections are interpolated to grid k, 
and to which residuals from grid k a-re transferred. We 
call £ "the predecessor of k", and k "a successor of 
£". Each subgrid, except for the coarsest, has exactly one 
subgrid defined as its predecessor, and may have any number 
of successors. The domain on which each subgrid is defined 
is always contained in, or coincide with, the domain of its 
predecessor. Thus, the set of subgrids is arranged 
logically in a tree, instead of the linear ordering we had 
before. 

All these subgrids are still uniform, and can still 
easily be handled (created, extended, displayed, etc.) by 
the system mentioned above. Except that some of the inter- 
polation routines required for this more general situation 


- 46 - 


rrm m . nni 

^ig. 3. A piece of non-uni form , boundarg-lager type grid 
(A) and the uniform rectangular subgrids' it is made of 
(B/ Cf 13 f E) • 
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are still missing in ^^UGTAPE (1978a); they have been 
prepared for MUGTAPE (1978b). 

In case the thin layer is not along coordinate lines, 
some local coordinate transformation is required. This 
technique is explained in Section 3.5 below. 

3.3. Difference Equations and Mu2ti~Crid Processing 

The difference equations and their multi-grid solution 

for the composite structure are explained in Section 2.2 of 

Brandt (1977b). The main idea is that a f ine-to-coarse 

H H 

correction function correcting the G equations, 

(see (2.16)-(2.17) and the subsequent discussion) can be 

computed wherever the finer grid G^ exists, or, more pre- 

cisely, at any interior point of G which is also an 

interior point of its successor G . At all other points 

the original coarse grid equation (2.3) can be used. From 

this point of view it is clear that we can have various 

patches of finer subgrids ("successors") thrown over various 

' H 

desired subdomains of G ; the finer-subgrid accuracy will 

be established in the equations of such a subdomain via the 
H . 

correction. On any part of such a subdomain a patch of 
still-finer subgrids can be defined, etc. 

The multi-grid process proceeds essentially as before. 

If, for example, we regard the coarsest subgrid as level 1, 
and all the successors of level j as level j + 1, we 
could use exactly the same algorithm; e.g., the one shown in 
Table 1 above. Except that now, each operation (relaxation, 
finc-to-coarsc residual transfer, etc.) at each level is 
performed not on one subgrid, but on the sequence of all 
subgrids of that level in some preassigned order. Important 
improvement: Relaxing each cycle progressively smaller 

parts of the coarser grids. 
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Note that successors of a given subgrid may geomet- 

rically overlap. All is needed in such a case is to set 

priority relations between the successors, to tell which 
H H 

correction 't. applies at those G poincs which are 
h 

interior points of more than one successor. Such priority 
relations are automatically implied by the ordering of sub- 
grids within each level. 

An important advantage of this structure is its flexibil- 
ity. One can add more and more patches of Increasingly finer 
subgrids, where and when they are needed. One can also dis- 
card some such patches. Notice, however, that even when a 

piece of finer grid, G^ say, is discarded one can still 

H H 

retain its correction t. in the G equations. (This 

El 

leads to multi-grid procedures which require only a small 
memory. See Section 2.4 in Brandt (1977b).) 

Another important advantage of the outlined structure is 
that our difference equations are defined on uniform grids 
only (patched together by the usual multi-grid interpola- 
tions). Such difference equations on equidistant points 
are simple and can be read from small standard tables, while 
on general grids their weights would have to be recomputed 
(or stored) separately for each point, entailing very 
lengthy calculations especially for high-order approxima- 
tions. Thus, our system facilitates the use of high and 
variable (adaptive) order of approximation. 

Still another advantage is that relaxation sweeps, too, 
are done on uniform grids only. This simplifies the sweep- 
ing and is particularly essential where s>'mmetric and/or 
alternating-direction line relaxations are required for 
obtaining high smoothing rates. 
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3.4. Remark on High-Order Approximations 

An efficient and convenient way of using high-order 
difference equations, especially when the order is adaptible 
(see Section 3.6), is by the well-known technique ^suggested 
by L. Fox) of "deferred correction" (see, e.g., Pereyra 
(1968), Lentini and Pereyra (1975) and Stetter (1978)). 
Simply, before starting a multi-grid cycle for improving an 
approximation u , add to the right-hand side of the fine- 
grid equations (2.2) the correction 


h f h. , b h . h. « h h . h. / h »h. 

Op(x ) - L u (x ) - LpU (x ), (x c G ) , 


(3.3) 


where is the higher-order (order p) operator. Then 

proceed with the multi-grid cycle as usual. The roll of 

is similar to the roll of the f ine-to-coarse correction 
H h 

r. . We can thus call c the high-order-to-low-order 
h p ^ 

("deferred" , or "defect”) correction. 

A certain amount of work is saved If p is advanced in 
steps; e.g. each multi-grid cycle advance p by 1 or 2. 
In the adaptive procedures described below, p is always 
advanced gradually. 

Note that, since the multi-grid cycle operates with the 
original operator L , no new routines (such as relaxation 
routine) should be added, and the same multi-grid efficiency 
is obtained as in solving lov?-order equations. Except that 

the number of multi-grid cycles may increase linearly with 

, 2hi h _ „p 

p, since t /t *= 2 . 


3.5. Local Coordinate Transformations 

Another dimension of flexibility and versatility can be 
added to the above system by allow’ing each subgrid to be 
defined in terms of its own set of coordinates. 
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Near a boundary or an interface, for example, the most 
effective local discretizations arc made in terms of local 
coordinates in which the boundary (or interface) is a coor** 
dinate line. In such coordinates it is ea-y to formulate 
high-order approximations near the boundary; or to introduce 
mesh-sizes which are much smaller across than along the 
boundary layer (see Section 3.2); etc. 

Usually it is easy to define suitable local coordinates 
(see below) , and uniformly discretize them, but it is more 
difficult to patch together all these local transformations, 
especially in an ad>‘)ptlble way. In the above structure, 
however, this difficulty does not arise, since we can 
Introduce independent and overlapping patches of "successor" 
grids . 

Each set of coordinates will generally have more than one 
subgrid defined on It, so that (1) local refinement, in the 
style of Figure 2 and/or Figure 3 above, can be made within 
each set of coordinates; and (ii) the multi-grid processing 
retains its full efficiency by keeping the mesh-size ratio 
between any subgrid and its predecessor properly bounded 

Since local refinement can be made within each set of 
coordinates, the only purpose of the coordinate transforma- 
tion is to provide a certain grid direction, i.e., to have 
a given manifold (e.g., a piece of boundary) coincide with 
a grid hyperplane. We can therefore limit ourselves to 
simple forms of transformations. For example, in 2-dimen- 
sional problems, let a curve (a boundary, an interface, etc.) 
be given in the general parametric form 

X » Xq(s), y*yjj(s), (0 1 s £ s^^) (3. A) 

%?here s is the arclength, i.e.. 
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( 3 . 5 ) 



i 



I 


To get a coordinate system (r,s) in which such a curve 
vJ.ll be a grid line, we can always use the transiormation 
x(r,x) - Xq(s) - ry^(s), y(r,s) - y^Cs) + rx^(s) . (3.6) 


Hear the given curve (r ■ 0) this transformation (a 
special case of transfonaations discussed in Starlus 
(1977a)) is orthogonal, owing to (3.5), and transforms any 
small h x h square to another h x h square. 

The main advantage of this transformation is that it is 
fully characterized by the single-variable functions Xq(s) , 
yQ(s). These functions (together with >:q(s), yQ(s) and 
q(s) ■ x^/y^ * “y^/xQ) can be stored as one-dimensional 
arrays, in terms of which efficient interpolation routines 
from (x,y) grids to (r,s) grids, and vice versa, can be 
programmed once for all. (Such a general routine, however, 
is not easy to program, and is still missing in MUGTAPE 
(1978).) The difference equations in (r,s) coordinates 
are also simple enough in terms of these arrays. For 
example, by (3.5-6), 




yo 


TZ • (3.7) 


3x "^0 9r 1 + rq 9s’ 9y 0 9r 1 + rq 9s 

Hence we can easily approximate the (x,y) derivatives by 
(r,s) finite-differences, with numerical values of Xq(s), 
>’q(s) and q(s^ directly read from their stored tables. 

(No interpolation is needed if the tables contain values 
for s points which correspond to grid lines and half-way 
between grid lines.) 

Such a system offers much flexibility. Precise treatment 
of boundaries and interfaces by the global coordinates is 
not required, since along boundaries the global grids arc 
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only correction to the local ones. The local coordi- 

nates are easily changeable (changing only the one<- 
dlmenslonal tables of ^Q*TQ>^Q>yQ*^) can therefore 

be adapted to a moving interface. 

The main difference between this structure and the one 
used by Starius (1977a), (1977b) and (1978) is that the 
boundary grids are completely embedded in the global grids 
(their predecessors), allowing a fast multi-grid solution 
of the equations. Also, since we have the multi-grid method 
for local refinement, the coordinate transformation is used 
only for orienting the grid, hence only the simpler trans- 
formation (3.6) is needed, allowing simpler differencing 
and interpolations. 

Another variant of this procedure is required in case 
Che location of the curve (interface, shock, etc.) is not 
fully defined. For example, a solution may include many 
shocks, some weaker and some stronger, and it is hopeless 
to try to recognize where a shock occurs, let along deter- 
silne its exact curve. The usual procedure is to let the 
shocks develop by themselves, e.g., by adding some avtificial 
viscosity which spread shocks over several mesh-sizes. 
Sometimes, however, this procedure is unacceptable because 
too much artificial viscosity is used near strong shock (and 
because of other reasons) . We like to have a procedure 
which will automatically use smaller mesh-sizes near stronger 
shocks. This will be done by the general adaptive procedure 
(Section 3.6 below) if we choose the error functional E 
so that it contains some measure of the artificial viscosity. 
In order to obtain full efficiency, however, we like the 
procedure to be able to produce mesh-sizes which arc much 
smaller in one direction (the direction perpendicular to the 
shock) than in the other. We therefore need a structure for 
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adapting Che grid orientation, too. Notice that on some 
coarse grids the orientation is imnaterial; the finer the 
grid the more precisely its orientation should be chosen. 
Hence, the direction can be refined succccsively , together 
with the mesh-sizes. An example is drawn and explained in 
Figure A. We sec that in this method the more general 
transformation (3.6) is not needed; only rotations are used. 
Hence the difference equations may be as simple as in the 
original (e.g., cartesian) coordinates. This method may 
therefore be preferable even in cases the curve (e.g., 
boundary) is known. 


3.6. AdtptMtion Technigves 


The flexible organization and solution process, described 
above, facilitate the implementation of variable mesh-size 
h(x) and the emplo>Taent of high and variable approximation 
order p(x). How, then, ?ire mesh-sizes and approximation- 
orders to be chosen? Should boundary layers, for example, 
be re.solved by the grid? What is their proper resolution? 
Should high-order approximation be used at such layers? How 
does one detect such layers automatically? In this section 
we survey (for more details, see Brandt (1977a), Chapters 8 
and 9) a general framework for automatic selection of h(x), 
p(x) and other discretization parameters in a (nearly) 
optimal way. This system automatically resolves or avoids 
from resolving thin layers, depending on the goal of the 
computations, which can be stated through a simple function. 

As our directive for sensible discretization we consider 
the problem of minimizing a certain error estimator E 
subject to a given amount of solution work W (or minimiz- 
ing W for a given E. Actually, the control quantity will 
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Figs. 4 a and 4B. Grid orientation around an interior thin 
lager. The two coarsest levels (A) have the usual orienta~ 
tion 0. The next level (B) has 3 orientations: 0, ^ and 

(the later is not applied here). The next level (not 
shown) would have 7 orientations: 0, etc. The 

successors (refinements) of a grid will always have either 
the same orientation or one of the two closest ones (e.g. , 
each successor of the --oriented grid in B will have 

, ir IT 3 tt, 

oraentataon — or -r-J . 

4 9 o 
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be neither E nor W, but their rate of exchange) . This 
optimization problem should of course be taken quite loosely, 
since full optimization vould require too much control work 
and would thus defeat its own puirpose. 

The error estimator E has generally the form 

E « / G(x)i(x)dx . (3.8) 

U 

where r(x) = r(x,h,p) is the magnitude of the local trunca- 
tion error at x (see (2.18)). G(x) 0 is the error- 
weighting function. It should in principle be imposed by 
the user, thus defining his goal in solving the problem. In 
practice G(x) serves as a convenient control. It is only 
the relative orders of magnitude of G(x) at different 
points X that really matter, and therefore it can be chosen 
by some simple rules. For example, if it is desired to' 
compute £-ordcr derivatives of the solution up to the 
boundary then 

G(x) , (3.9) 

where d^ is the distance of x from the boundary, and m 
is the order of the differential equation. 

The work functional W is roughly given by 


W *= / dx , (3.10) 

0 h(x)° 

wiiere d is the dimension and h is therefore the number 
of grid points per unit volume, w «= w(p) is the solution 
work per grid-point. In multi-grid processing, for p ^ 6 
the work depends mainly on the number of cycles (cf. Section 

3 

3. A), hence w » w^p. For (unusually) large p, v *= 0(p ) 
since evaluating at each cycle involves 0(p) 

0(p) arithmetic precision. 


terms and 


Treating h(x) and p(x) as continuous variables, the 
Euler equations of minimizing (3.8) and (3.10) can be 
written as 

h 

C ll + 2 e1|eI . 0 (3.11b) 

3p 

where X is a constant (the Lagrange multiplier), repre- 
senting the marginal rate of exchanging optimal accuracy 
for work: X = -dE/dW. The * sign in (3.11b) should be 

replaced by ^ at points x where p attains its minimal 
allowable value p^ (usually 1 or 2) . In case we use 
fixed— order difference equations (adapting h(x) only, p 
is fixed in advance), equation (3.11b) should be omitted. 

If constant order is used (p is constant over the domain, 
but instead of being fixed in advance this constant is to 
be optimized) equation (3.11b) is replaced by 

|£+x|H-o. (3.12) 

3p 3p 

In principle, once X is specified, equations (3.11) 
determine, for each x c fi, the local optimal values of 
h(x) and p(x), provided the truncation function 
t(x,h(x) ,p(x)) is fully known. In some problems the 
general behavior of t(x,h,p) near singularities or in 
singular layers is known in advance by some asjTnptotic 
analysis, so that approximate formulae for h(x) and p(x) 
can apriori be derived from (3.11). More generally, how- 
ever, equations (3.11) are coupled with, and should there- 
fore be solved together with, the given differential equa- 
tions. Except that (3.11) is solved to a cruder approxima- 
tion. This is done in the following way; 
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In the multi-grid solution process (possibly incorporat- 
ing a continuation process) , incidentially to the stage of 
computing f from u , u and f (see Figure 1 above, 

2h corresponding to level k and h to k + 1) we can 
get an estimate of the decrease in the error estimator E 
introduced by the present discretization parameters. For 
example, the quantity 


-AE(x) = G(x) lT(x,2h,p) - T(x,h,p)I (3.13) 

« G(x)1t^^| 

= G(x)lf2h(x) - F^^(x) 1 

(cf. (2. 17) -(2. 19) above) may serve as a local estimate for 
the decrease in E per unit volume owing to the refinement 
from 2h to h (cf. (3.8)). Each such decrease in E is 
related to some additional work V (per unit volume) .. For 
example, that refinement from 2h to h required the 
additional work (per unit volume; cf. (3.10)) 


h** (2h)‘^ 


(3.U) 


Ve say that the present parameter (h in the example) is 
highly profitable if the local rate of exchanging accuracy 
for work Q = -AE/AW is much larger than the control para- 
meter X. 

More sophisticated tests may be based on assuming t to 
have some form of dependence on h and p. Instead of 
calculating Q for the previous change (from 2h to h in 
the example) we can then extrapolate and estimate the rate 
Q for the next change (from h to h/2), which is the 
more appropriate rate in testing whether to make that next 
change. The "extrapolated test" is not, however, much 
different in practice, and may actually be equivalent to 
testing the former Q against another constant X. 
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1q deciding whether and where to change the discretiza- 
tion, ve adopt rules which stablize the adapt.' ve process. 

For example, a change (e.g., refinement from n to h/2, 
which in practice (see Sec. 2.2) means an extension of the 
uniform grid with mesh-size h/2) is introduced only if 
there is a point vrhere the change is "overdue" (e.g., a 
point where Q > 15X) . But, together with each point where 
the change is introduced, it is also introduced at all 
neighboring points where the change is just "due" (e.g., 
where Q > 3X) . 

We can use the Q vs X test to decide on all kinds of 
other possible changes, such as: Changing the order p 

to p + 1 (or p - 1) ; or changing the computational 
boundaries (when the physical domain is unbounded) ; or we 
can use such a test to decide whether to discard some terms 
from the difference operator (such as the highest order 
terms in some regions of a singular-perturbation problem) ; 
or to decide on unisotropic changes in h and p (e.g., 
changing Ax to Ax/2 without changing Ay); etc. 

In case of optimizing a global quantity (such as constant 
approximation order; cf. (3.12)) similar tests can still be 
applied, except that AE and AW should of course be 
measured globally (summing over the entire region) instead 
of locally. 

The computer work invested in the tests is negligible 

compared with the solution work itself. The measure (3.13) 

2h h 

is taken only on G g , and the stage of computing it 
occur only once per several relaxation sweeps on G . 
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A. multi-level adaptive solutions to singular-perturbation 

PROBLEMS: CASE STUDIES 

To get a transparent view of the discretization patterns 
and the accuracy-work relations typical to the adaptive 
procedures proposed above, when applied to singular perturba- 
tion problems, we consider now several cases which are 
simple enough to be fully analyzed. The simplicity of the 
solution, It should be emphasized, is not used in the solu- 
tion process itself. 


4.1. Optimal Discretization of One-Dimensional Case 


Consider the 2-point boundary-value problem 
0 in 0 < X < 1 , 


dT . dU 
dx 


(4.1) 


with constant c > 0 and with boundary conditions U(0) 

— x/s 

and U(l) such that the solution is U *= e .An elliptic 
(stable) difference approximation to such an equation can 
be central for c _> 2h but should be properly directed for 
small e. (See Section 5 below.) In either case, the 
truncation error behaves like 


x(x,h,p) 

where y is 
depends on 
dependence .) 
G(x) H 1 


lyej 


p g-x/c 


(4.2) 


a constant close to 2. (Actually, y slightly 
p; see (5.14). For simplicity, we neglect this 
For the error weighting function we choose 
, (4.3) 


which would be the choice (see (3.9)) when one is interested 
in computing boundary first-order derivatives (correspond- 
ing, e.g., to boundary pressure or drag, in some physical 
models). Then, assuming w(p) w^p, the optimization 
equation (3.11a) yields 
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and (3.11b) therefore becomes 


if P 


h > ■^ c if p 


e - * ^0* ‘ - e " **0 

From (4.2), (4.4) and (4.5) it follows that 

h ■ c, p ■ log 


, - 1 for 0 < X < Xrt . 

AWq c “0 


(x“X.)/(ep.+E) 

ai X - - ^ 


c e 


for Xq<.x < 1 » 


where 


*0 ■ • Po • » • 

If Xq 1 then (4.6) applies throughout, hence 

1 w.p w.e , 

W ■ f — ^ dx ^ ^ ^ ^ 

I ^ 


(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 


YE AWq " ^ " 2e ^ * 


1 Xw e 

E “ / Tdx - — — 
n Ye 


7 exp(- W - . 


V G 
0 


(4.9) 


(4.10) 


and the condition 2. i itself becomes, by (4.8)-(4.9), 


1 


(4.11) 


0 2e 

Thus, if W satisfies (4.11), E converges like (4.10). 
That is, for large values of c, the total error E 
decreases exponentially as a function of the overall work V. 

Notice that when t is large no boundary-layer is 
formed, and the mesh is uniform. Note also that the optimal 
mesh-size h * Ye/e is independent of the work W (or the 
exchange rate X). That is, when more work can be afforded, 
it should not be invested in refining the grid, but in 
increasing the approximation-order p. In the MLAT 
processes described above p will automatically increase 
by 2 (or by 1, if non-central approximations arc used) 
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every multi-grid cycle, until the desired accuracy (or the 
work limit, or the prescribed exchange rate X) is reached. 

4.2. Boundari'-Layer Resolution 


For small c, however, the mesh-size h * c is 
impractical. Indeed, in the optimal discretization (A.6)- 
(A.7), for small t we get small x^, and an "external 
region" x^ ^ x < 1 is formed where the mesh-size grows 
exponentially. The small mesh-size is used only to resolve 
the boundary layer. In this simplified problem the solution 
away from the boundary layer (i.e., for x » c) is practi- 
cally constant, so that indefinitely large h is suitable. 
Usually h will grow exponentially, as in (4.7), from 
h ■ ye/e to some definite value suitable for the external 
region (the optimal mesh-size of the reduced problem) . In 
the transition region we have p ■ p^, i.e., the minimal 
order of differencing is used in the region where h 
changes. This may be useful in practical implementations. 

From (4.6)-(4.8) and (4.2)-(4.4) we get, for small c. 


W 


1 w-p w_e 

[(log ~ 
h 2y ^ Xw. 


1)^ + 2P(, + 


(4.12) 


1 Xw e 

E » J tdx = — — log ■— , (4.13) 

0 ^ 

where exp(-l/t) and similar terms are neglected. Using 
(4.12) we can express X in terms of W and substitute 
that expression is (4.13). In reasonable calculations 
W Wq, and then the relation simplifies to 


Thus, essentially - 



J J 


(4.14) 
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For small values of c, the total error F. decreases 

1/2 

exponentially as a function of V , where this rate is 
independent of c and does not deteriorate as e 0. In 
principle, th.'.s rate is better than O(h^), for any fixed p. 

Notice that here h docs depend on W (or X), but 
only in some transitional layer. In the inner part of the 
boundary layer h « ye/e still holds, while away from that 
layer h tends to the optimal mesh-size of the reduced 
problem. (If the reduced problem is itself regular, its 
optimal mesh-size will be determined by the "local scale" 
of that problem. This scale is independent of V, as it 
is for example in Section A.l above for the case of moderate 
c. That scale is too small to resolve only when the reduced 
problem is singular) . ffhat depends on the total computa- 
tional work is the distance from the wall at which 

the meshsize starts to grow exponentially . In fact, from 
(4.8) and (4.12) we see that 
, 1/2 

. (4.15) 


JOL y 

"o® • 


Defining the computational boundary layer as the region 
where h < h^ for some h^ independent of c, the width 
'■’CBL layer is, by (4.7), 

''CBL 7 • (4.16) 

Another, quite obvious but interesting observation can be 
made at this junction, based on (4.11) above. Even for 
small e, if W is sufficiently large then the exponential 
relation (4.10) holds. Hence, in an asjTnptotic theory for 
W ^ 00 (corresponding to asymptotic theory for h -► 0, 
which is so common in numerical analysis) the relations of 
this section, which undoubtedly dominate the numerical 
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process at reasonable values of M, would not be sceni 
What we should be interested in are values of W which are 
large but independent of C/ and therefore not large 
compared with negative powers of c. 


4.3, Fixed-Order and Constant-Order Discretizations 

The optimal approximation-order p calculated above 
varies with the location x. This is not essential. 

Indeed, if p is fixed then (3.11b) is omitted, but (3.11a) 
and (4. 2)- (A. 3) still imply x » Xw(p)/(hp), and hence 


[Xw 
^ YP 


x/(cp+t) 


ee 


Hence, for small c, 




1 

p+1 


Yw 


1 w 
P i h p 


Xp + 1) 


f 


yV 


i(p + l)wj 


-P 


(4.17) 


(4.18) 


(A. 19) 


so that the 
E - O(h^) 


Thus, E = CW”P with C independent of e, 
convergence order is p (analogous to error 
when h is constant). Tha variable mesh-size (4.17) keeps 
the convergence rate essentially unimpaired by the singfular 
perturbation, even though convergence is considered of the 


first derivative up to the boundary. 

The constant approximation-order p need not of course 
be fixed in advance. It may be optimized just as well, 
using global tests as mentioned above. From (4.19) it 
follows that the minimal E for a given W is obtained 
when p satisfies 
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(P + l)w(p) 


(4.20) 


1 + 


pw*(p) , 
w(p) 


» 


and for w 
E 




cht: total error will be 


, „.l/(£+l) 
(yW) "exp 


1 + £ 

f \ 

l/(£+l)‘ 

e 

■''o- 

< 


(4.21) 


For I ■ 1 thifi rate Is alnost the same as (4.14). Thus, we 
do not lose much by using constant, optimized p, which, 
on the other hand, may be considerably simpler to program. 

From (4.17)“(4.18) and (4.20) we see that the width of 
the coD^utational boundary layer is now 

''CBL ** c ’ (4.22) 

For small c this is (p + 1 )/(Pq + 1) times wider than 
the variable-order case (4.16). 


4.4. A Case of Skipping the Boundary-Layer 

To see the effect of choosing different error-weighting 
functions, consider again problem (4.1), but vd.th the choice 

C(x) = X . (4.23) 

This will be the choice in case one is interested in 
approximating U only, not its derivative , and to approxi- 
mate it in the sense. By substituting (4.2) and (4.23) 

into (3.11), and solving for p, we would obtain 

P “ log ^ f (4.24) 

- 2 • 

0 

For bounded (independent of e) X and sufficiently small 
c, this p is smaller than p^. Hence P “ Pq should 
replace (3.11b) and substituting (4.2) into (3.11a) we 
actually get 
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Thus, for bounded X and sufficiently small c, h >> c 
everywhere, so that the boundary-layer is not resolved by 
the grid. 

In fact, since h >> c, all interior grid points lie in 
a region where the rate of convergence would normally be 
determined by the reduced equation (see Section 4.5). In 
our simple ex.'unple (4.1), the reduced equation has the 
trivial solution U = 0, and accordingly E 0 as c 0, 
for any fixed W (or X). This can be verified from (3.8), 
(4.2) and (4.25). 

In case one is interested in approximating U in the 
sense, a precise choice of the error-weighting function 

is 

G(x) • 1 - e . (4.26) 

With this function, solving (3.11) for p we would get 


P(x) 



(4.27) 



max p(x) - p(c log 2) - log 1 , 

and hence, for X reasonably small, p(c log 2) > p^. 

y 

Therefore, around x = e log 2, we again have h “ -^e. 

■ Beyond this point (for x > c log 2) the discretization 
pattern is essentially the sane as in Section 4.2 above 
(since G(x) is essentially the sane). Before this point 
(x _< c log 2) we have h(x) > x, so that in practice we 
do not have there more grid points. Thus, the grid through- 
out is essentially as in Section 4.2. Similarly, for 
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J 


constant p the mcsh-size distribution will be as in 
Section A. 3. The accuracy-work relation, too, is essentially 
as before. 


4.5. Remarks on More General Problems 


For pe-.ieral problems it is of course impossible to find 
apriori the relation between work and accuracy tlat would 
result from multl~level adaptive solution processes. In 
fact, in most non-trlvial cases, an optimal (or nearly 
optimal) choice of discretization parameters (h(x), p, 
etc.) is not known in advance, since it depends on the 
particular solution. This is exactly why adaptive tech- 
niques are needed. Nevertheless, in this section we will 
try to indicate that the simple relations described in 
Section 4. 2-4. 4 are typical to manv, perhaps most, singular- 
perturbation problems, even in complicated, high-dimensional 


problems . 

Consider first a more general one-dimensional, 

coefficient equation of the fora 

m-n^(m) ^ ^ ^m-n-Iy(m-l) ^ ... ^ „(n) 

m-1 n 


+ a + ••• + U 

n-1 


( 0 ) 


constant- 

(4.28) 


normalized so that 

VM . (‘-29) 
is a solution. And assume the boundary conditions are such 
that (4.29) is actually the solution. The truncation error 
is then approximately (see (5.14)) 



(4.30) 
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where y is again a constant (slightly different than 
in (4.2), but still Y 2 for larger p. We again 
neglect the changes in y). 

If we are interested in computing near x = 0, 

the error-weighting function for small x should behave 
like 


G(x) . (4.31) 

For the adaptive process, the multiplicative constant in G 
is immaterial. For our convergence estimates, however, the 
correct order of e should be used. The behavior (4.31) 
results from the observation that if ^'^^^(x) = 6^(x) 
in (0,1) and V(0) = V(l) = 0, then, for 0 < x < C it 
tx) = 0(e ^) . An additional factor appears 

in (4.31) since we are interested in measuring relative 
errors in u^'^^(x), and by (4.29), u^^^(x) = 0(c ^) for 

X < 0(c). 

For j = m - 1, Gt is the same as in the special case 
discussed in Sections 4.2 and 4.3 above. Therefore exactly 
the same discretization parameters and the same accuracy- 
v?ork relations will followr. For smaller j, the accuracy- 
work relation cannot get worse; it may even improve, 
depending on the norm used (cf. Section 4.4). 

In more general singular-perturbation problems, the solu- 
tion U(x) can be written as a sum of a function U(x) 
which tends uniformly to the solution Uq of the reduced 
problem, and bouedarj'-layer terms, each of which behaves 
like (4.29) above for some suitable c. (See O'Malley 
(1974).) Consider the case of a fixed order p, as in 
Section 4.3 above. Let h^(x) be the mesh-size distribu- 
tion optimal (at some given X) for the reduced problem, 
and h^(x) the optimal distribution in case the solution 
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contain only Llie i-th boundary tern (i = l,2,...,o). Let 

and denote the corresponding total error and overall 

work (i = 0,1,2 o) . For any given solution (containing 

all terms) choose 

h(x) = min h.(x) . (4.32) 

O^ij^a ^ 

Then clearly E ^ F and W ^ I W^. In an optimal choice 
of h, E vill be even smaller (for the same value of W; or 
vice versa). Hence, essentially. 

The convergence rate behaves either like one of those 
described above for the boundary terms (Section 4.3), or like 
the convergence rate of the reduced problem , whichever is 
si ower . 

The situation is a little more complicated for variable 
p, but we saw before that we don't ]oose much by using a 
constant p. Moreover, the optimal p (4.20) does not 
depend on c and can therefore serve uniformly for all the 
boundary-layer terras. 

Similar convergence rates should be expected in higher'- 
dimensional problems, too. To see this, examine the 
behavior near some portion of the boundary. Assuming our 
computations use boundary coordinates (see Remark below) , 
we can regard the boundary as = 0. Using Fourier 
transformation in all but the x^ coordinate, the solution 
u can again be written, in many cases, as a sum of u 
(tending asymptotically to the reduced solution U^) and 
boundary-layer terms behaving like (4.29). Assuming also 
that the only significant adaptation of mesh-size is needed 
in the x direction (i.e., perpendicular to the boundary), 
we may repeat the above argument, using (4.32), and arrive 
at the same conclusion. 
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Remark about boundary coordinates : "Boundary Coordinates" 

is a coordinate system in which the boundary is contained, 
at least iocally, in a coordinate hyperplanc (c.g., 

{x^ = 0}). Ir. Section 3.5 above it is explained how to 
construct and use such a coordinate system in multi-grid 
processes. For the finite-difference equations it is 
important to use a grid along such boundary coordinates. 
Otherwise it is impossible to simultaneously use small mesh- 
sizes in the direction perpendicular to the boundary and 
large ones in the other direction(s) , as required for 
obtaining efficiencies similar to the one-dimensional ones. 

Thin transition layers not on the boundary, such as turn- 
ing points in ordinary differential equations or contact- 
discontinuities and shocks in higher dimensions, are lirtely 
to be treated by multi-level adaptive techniques as 
efficiently as the boundary-layer cases analyzed above, 
since the procedures did not assume any apriori knovcledge 
concerning the location of the layer. The layer is dis- 
covered, and if necessary resolved, by the numerical process, 
using general and automatic criteria. The only difficulty 
is, in higher dimensional problems, to get a coordinate 
system in which the interne 1 layer is a coordinate hyper- 
plane. To a suitable approximation, however, this can be 
done, using the procedure described in Figure 4 above. 

Not all singular-perturbation problems can efficiently be 
solved by the above techniques, of course. For example: 
problems with highly oscillatory solutions, such as the 
Helmholtz equation 

e^Au + u = 0 . (4.33) 

In usual norms, this problem is not uniformly well-posed . 
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That is, the change in the solution caused by a certain 
change in the data is not uniformly bounded: it may increase 

indefinitely as c 0. Such problems should be reformulated, 
using other variables and norms, so as to make them unifoniily 
well-posed. (See Section 5 in Brandt (1978b).) 

5. UNIFORMLY WELL-POSED HIGH-ORDER DIFFERENCE EQUATIONS 

An extended version of this section appears as Section 5 
in Brandt (1978b) . It discusses the concepts of well- 
posedness and uniform well-posedness, ellipticity and uniform 
ellipticity, and their significance for singular-perturbation 
problems in general, and for their multi-level solutions in 
particular. Closely related are the extensive theoretical 
investigations of Frank (1978 and references therein). 

Related preliminary’ observations were made in Brandt (1976) . 

Here we summarize some of the more practical aspects. 

5.1. General Remarks 

In approximating potentially singular-perturbation equa- 
tions it is essential to ensure that the discrete problem 
is uniformly well-posed (uniformly stable) not only with 
respect to the mesh-size (h) , but also with respect to 
the singular-perturbation size (c) . That is, in suitable 
norms, a small change in data should cause a small change 
in the solution, uniformly in both h and e. For this 
to be possible, the original differential problem should be 
uniformly stable (in e) . This, however, is not sufficient. 
Innocent-looking difference approximations L may easily 
be uniformly stable in h (i.e., for any fixed e) , but 
not jointly in h and c. In such cases satisfactory 
approximations will still be obtained by sufficiently small 
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h, but that h will have to be small compared with e (or 
some power of c), and hence too small to be practical. In 
particular such mesh-sizes are unacceptable (even for 
moderate c) for the coarse grids of a multi-level sturcturc. 

There are no general procedures to construct uniformly 
stable difference approximation; nor even general procedures 
to check uniform stability of given difference schemes. 

This is in fact already true for the differential equations. 
But there are some important classes of uniformly stable 
operators and some practical ways of constructions. 

For various boundary-value problems to be well posed it 
is required that the partial-differential operator (2.1a) 
is elliptic, i.e., that the homogeneous syste.*n of equations 
has no non-constant periodic solution. This, together with 
appropriate boundary conditions, ensures w’ell-posedness. 
Similarly, the difference operator (2.2) can be defined as 
elliptic if there is no periodic such that L^U = 0. 

For scalar operators (q = 1) such a definition was 
introduced by Thomee (1964), and various results related to 
the stability of such operators were proved by him and by 
Thomee and Westergren (1968) . Many more results were 
published in conjunction w^ith finite-element formulations, 
which yield scalar or vectorial elliptic difference operators . 
(See Ciartet (1978).). Most of these results, however, 
hold only for sufficiently small mesh-sizes, and are there- 
fore not directly applicable in the present context (where 
"sufficient-small" means smaller than e) . Slightly different 
notions of ellipticity are needed. The most useful for 
applications is, perhaps, the following. 

R~Ellipticity. Assume the q q difference operator 
in (2.2) has constant coefficients, and let 
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h = (h, d ,) , where h. :^s the mesh-size in the x. 

1 B j J 

coordinate. Denote 


6 = (0j^}...,6^), 

je 1 = max( 1 


6-x/h 




• • • + 6 



(5.1) 

(5.2) 


Then, for any constant q-vcctor V 

j^h^ie-x/hy . , (5.3) 

w’liere A is a q q matrix, easily obtained by replacing 

in the matrix each h. -translation with the complex 

^ h 

function exp(iO.). B is called the matrix-symbol of L . 

^ h 

The difference operator L is called R-elliptic if 

fle v'^B(e,h)V >0 for all 0 < le] £ and all (5.4) 


real q-vectors V # 0 


An operator with variable coefficients is called R- 

elliptic if the frozen-coefficients operator at every point 
is R-elliptic. A nonlinear difference operator i» called 
R-elliptic if the corresponding linearized operator is R- 
elliptic (which may depend on the solution around xdiich the 
linearization is taken) . 

This is not a complete definition of ellipticity. For 
example, if ar.d are R-elliptic, then 

not necessarily also R-elliptic. But the definition gives, 
on one hand, a concept much more general than the special 
case of positive-type operators (w’hich trivially satisfy 
(5.4)); in fact, a definition general enough for almost all 
scalar equations. On the other hand, the definition has some 
nice properties. One property is that it restricts the 
location of the operator, while Thomee's definition allows 
any translation to be added to the operator (which of course 
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cannot be permitted in discussing finite mesh-sizes, because 
for example, it allows the tvro difference equations at two 
neighboring points to coincide, i.e., to be just one equa- 
tion). This restriction is essential in discussing relaxa- 
tion schemes, where a relation is required between each 
difference equation and the point at which it is relaxed. 
Another nice property is that the sum of R-elliptic operators 
is clearly also R-elliptic. One can therefore construct 
R-elliptic operators one term at a time. 

We can use this property for singular -per turbation 
operators. If both the perturbed and the reduced equation 
are elliptic, the required uniform stability is obtained by 
constructing a difference approximation which is uniformly 
elliptic. A simple way to achieve this is to construct R- 
elliptic approximations to the various terms in the equation, 
so that R-elliptic approximation is obtained, in particular, 
for the reduced equation. 

5.2. Examples 

R-clliptic approximations, of arbitrary order, will be 
constructed in this section for the basic one-dimensional 
operators. Since this construction do not use any relation 
between terms, these approximations can be used as building 
blocks for approximating many ordinary and partial differen- 
tial operators. The approximations are constructed on 
uniform grids only. As shown in Section 3, this is all we 
need in a multi-grid environ.ment. 
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Using the operator notation 
fiu(x) « u(x+-|) - u(x--|), Vu(x) 
Au(x) “ u(x + h,/ - u(x), 
yu(x) « ~ [u(x+^) + u(x-j)], 


Du 


du 
dx * 


u(x) > u(x - h), 


(5.5) 


and the calculus of such operators (see, e.g., Dahlquist 
and Bjork (1974) p. 311) one can derive the expansions 

hH-vS I 3 (-tY (5.6) 

qcQ ’ 

(hD)^ = <5^ I (5.7) 

q=0 ^ 

where 

a 

®0 “ 2 hence (a^)^'^*^ ^ . (5.8) 

q 

From this ve find the following expressions for the 2s-order 
central approximations to the first and the second deriva- 
tives, and for the corresponding local truncation errors: 


1 

u'(x) = ^ p6 I 
q=0 

+ (-l)'®ag h^® 


a^(-{2)‘J u(x) 


^(2s+l) 




(5.9) 


, , s-1 a 

ru"(x) = (-6^) I ^ 


q + 1 


(-6^)'^ u(x) 


-(- 1 ) 
where C. 


a 

s s 


q=0 
2s (2s+2) 


(5.10) 


(^2) , 


1 

and ^2 3^® some intermediate points. 
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Let us now check these difference approximations for 

R~cllipticity. It is easy to see that the s^mibols corre- 

2 

sponding to p6 and to -6 arc, respectively, isinS and 

2(1 - COS0) . The latter is positive for all 0 < [©I £ tt. 

Hence all the above approximations to -u" (note the sign!) 

are R-olliptic. On the other hand any central approximation 

to either u' or -u' has purely imaginary symbol, and 

is therefore never R-elliptic. 

In various elliptic equations, -u" is added to au*, 

where a may have any sign. We therefore need to construct 

R-elliptic approximations to both u' and -u'. This is 

done by a'^ding to (5.9) an R-elliptic term of order high 

2s-l 

enough: To obtain an approximation 0(h ), add any 

12s 

positive multiple of the term — (-6 ) ; to retain the 

2s ^ 12s 

0(h^ ) order, add any positive multiple of -r V(-6 ) or 

1 2 s ” 

-•r A(~5 ) . These terms are R-elliptic since the symbol 

“ •*16 16 
of V and -A are 1 - e and 1 + e , respectively. 


< IT. 


The 


so that their real part is positive for 0 < |G. _ 

values of the positive multiples can be chosen so that the 
£ 

0(h ) approximation uses exactly £ + 1 points 
(£ = 2s - l,s). This gives the following R-elliptic 
approximations and truncation errors: 

s-1 . _ a 


± u'(x) = ^ pi I ) + 


+ (- 1 ) 
u' (x) 


a , 
s s-1 


q=0 
2s-l (2s) 


(-6^)®[u(x) 


(5.11) 


(43) 




s-1 


p6 


y a_(-6^)'l + 


q=0 


+ (-!)■ 


s-1 


+ a 


2h 

2s (2s+l) 


V(-6^)®lu(x) 


(5.12) 
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- (-!)■ 




h^s 


Observe that these operators are completely one-sided 

2 

(so called "upwind") only for the 0(h) and 0(h) 
approximations. One can describe the above formulae as the 
non-central operators closest to the central among all 
operators which use the minimal number of points. The one- 
sided operators using the same number of points are not R- 
elliptic (for orders higher than 2) . 

The error term. For theoretical purposes (as in Section 
4.1 above) it is more convenient to express the magnitude of 
the error terms in (5.10) and (5.13) in the forms 

2s 

(5.14) 


h 

2 

and 

f \ 

h 

1 



[yj 




respectively, and similarly for (5.11) and (5.12). It is 

^ and 2 


clear from (5.8) that both y*" 2 

s 


as s 

' 8 , 

Y^ does not 
s 


grows. In fact, as shown in Table 2, each 
change much with s, and for theoretical convenience we 
treat them as constants. 


TABLE 2 


s 

1 

2 

3 

4 

5 

-1 

a 

6 

30 

140 

630 

2772 

s 







1.22 

1.71 

1.86 

1.93 

1.97 


3.46 

3.08 

2.87 

2.73 

2.64 



High-ordcr approximations near boundaries laay pose a 
problem, since the above difference operators may need func- 
tion values at points which arc not on the grid. One way 
out is to impose this technical restriction on the adaptive 
process (Section 3.6) which will, as a result, choose to 
further refine toward the boundary. The refinement will be 
geometric, so that without using too many points (their 
number is proportional to the high approximation order 
desired in the interior), the mesh-size near the boundary 
will be small enough to allow low-order approximation. In 
this respect the boundary behaves like a singular curve. 
Incidentally, for certain error norms (correspondingly: for 

certain functions G) , lower order can be used (correspond- 
ingly: will be affected by the adaptive process) near the 

boundary without spoiling the global order of approximation. 
(Cf. Bramble and Hubbard (1962).) 

6. RELAXATIONS OTTH UNIFORM SMOOTHING RATES 

A full version of this section appears as Section 6 in 
Brandt (1978b). Here we summarize the main points through 
simple examples. 

The role of relaxation sweeps in multi-grid algorithms 
Is to smooth the error (Section 2.1). The efficiency of 
relaxation is therefore measured by its "smoothing factor" 
p and the corresponding "smooching rate" v « |log p| 

The smoothing factor is defined in terms of the local mode 
analysis. Namely, if p(6) is the convergence factor per 
relaxation sweep of the 6 Fourier-component (see for 
example (2.7)) then, 

y « max y(6), where isl «=. maxle.l . 


( 6 . 1 ) 



The range £ [oj £ ti, which Is somewhat arbitrary (see 
Section A. 3 in Brandt (1977a)), is chosen because these are 
the components which are too high to be seen on the coarser 
(2h) grid, so they cannot normally be reduced by the 
coarse-grid corrections. This definition seems to assume 
an infinite domain (where the Fourier expansion is made), 
but the behavior of such high-frequencies is practically 
independent of the domain. Thus, the smoothing rate v is, 
roughly, the number of relaxation sweeps required to reduce 
all high-frequency components by the factor 1/e. p and v 
can be calculated for any relaxation scheme by the MUCTAPE 
(1978a) routine SMORATE. A table of representative values 
is given in Brandt (1977a), pp. 351-352. 

For uniformly elliptic operators all (reasonable) relaxa- 
tion schemes have bounded smoothing rates. (See the general 
theorems in Section 3.1 of Brandt (1976).) For singular- 
perturbation problems, however, many relaxation schemes will 
j have V which grows indefinitely as the siae of the 

perturbation decreases (c 0) . That is, the convergence 
rates of some components 0 will not be bounded uniformly 
in c. One may sometimes still get a nice multi-grid 
process if those bad components have only a small contribu- 
tion to the error norms (sec Poling (1978)), but it is 
better and safer to use other relaxation schemes, with 
smoothing rates which are bounded uniformly in c. 

Two kinds of degeneracies w’ill usually occur in relaxing 
singular-perturbation problems. One kind occurs in the 
boundary layer, in dimension d > 2, when a highly stretched 
grid (as in Figure 3E) is used. To see the problem, 
consider the usual, pointwise Gauss-Seidel relaxation for 
the 5-point Laplace operator on a grid with aspect-ratio 
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a » h./h- << 1. Exaoininf, thr* component 0 * (0,—) for 

^ ^ - -2 ^ 
example, it is easy to see that v > .5a . Since a may 

be comparable to c , this smoothing rate may be extremely 

slow. The same slow rate will occur for any point-wise 

relaxation. If, in addition to the Laplace operator, the 

differential equation has also lower-order terms (as it does, 

of course, in singular-perturbation problems), the trouble 

still occurs, since on the finest grid the higlier-order term 

(the perturbation) is dominant. 

This kind of trouble can always be avoided by using 
Gauss-Scidcl lino relaxation. This means a relaxation in 
which we scan (cf. Section 2.1) not point by point, but 

line by line. At each line, all the values u^(x ) 
associated with that line are simultaneously replaced by new 
values which are computed so that they simultaneously satisfy 
all the difference equations associated wtith that line. In 
the above example the lines should be horizontal lines 
(Xj^ » const.), and the resulting smoothing rate will 
uniformly be v » 2/log 5, no matter how small a is. The 
same smoothing rate will be obtained generally in the 
boundary layer, provided line relaxation is employed with 
lines perpendicular to the boundary. 

Another type of degeneracy occurs on the coarser grids, 
where the reduced part of the equation dominates the smooth- 
ing process. The relaxation there should be one which is 
suitable for the reduced equation. Still more difficult nay 
be the case of intermediate grids, where both the reduced 
and the perturbation parts interact with the smoothing 
process. Consider for example the ordinary differential 
operator 


( 6 . 2 ) 


, ,, - , d^u du 
LU . + a ~ , 

dx 

8) 

and its lovcst-order ^ stable approxlcatlo: s (see Section 
5.2) 

= ~ {U^(x-h) - (2 + n)U^(x) + (1+ n)U*’(x + h)) (6.3a) 
h 

for a ^ 0 , 

lV = ~ {(1- n)U^(x> h) - (2-n)U*'(x) + U^(x + h)} (6.3b) 


for a £ 0 , 

where n * ah/c may be moderate or large,. Denote by y^(ri) 
and P_(n)i respectively, the smoothing factors for the 
forv3”d and backward Gauss-Scidcl relaxation of (6.3). 
Forward and backward refer to the marching direction, i.e., 
to the order in which the points x are relaxed. A 
straightforward calculation gives, for h ^ 0, 

1 + n I 


p^(n) “ P_(-n) 
y^(n) * y^(*-n) 


2 + n + il 

1 


(6.4) 


(6.5) 

which means that the 


|2 + n + i(l"*'ri) 

Observe that v.(n) -►I as n •♦ *, 

+ 

forward relaxation is not uniformly smoothing for a > 0, 

and should not be used. No "relaxation parameter" will 

help here. On the other hand in this case (a > 0) we have 
- 1/2 

P (h) <5 , so the backward relaxation has a very good 

uniform smoothing rate. The backward direction corresponds 
to the direction of convection , or the down stream direction, 
in physical problems modelled by (6.2). Generally, 
physical insight is an invaluable source for devising 
successful relaxation schemes. 


For a < 0 the situation is reversed: Backward relaxa- 

tion is useless at small c, but the forward relaxation 
has exccllenL (very small) smoothing rates. Slightly more 
difficult is the case where a = a(x) changes sign in the 
domain. In that case each relaxation direction will have 
slow smooth'* ng at some part of the dexnain. The good scheme 
then is symmetric relaxation^ i.e., sw’oeping forward and 
then backward. The smoothing factor (per single sweep) of 
this is 

y (n) = [y,(n)y (n)]'^^ = 

S T • 


1 + n 

3 + 3n + h + i(2+ri) 


j-/ ^ 


( 6 . 6 ) 


,- 1/2 


Hence _< 2/log 3, uniformly bounded for all values of 
n, positive or negative. Observe that, in fact, the 
larger is jnl the better is the smoothing rate. 

The same holds for singular-perturbation equations in 
higher dimensions: Very good smoothing rates are obtained 

by a proper choice of the relaxation marching direction. 

In some situations all marching directions should be employed 
successively if a uniform smoothing is to be achieved. 

This may require more sweeps per multi-grid cycle, which 
one w’ould like to avoid. We can, in fact, construct relaxa- 
tion schemes which have bounded smoothing rates even when 
marching against the local convection direction . These 
schemes necessarily belong to the following class. 

Distributed Relaxation . In classical relaxation v?e 
relate the unknowTi at a grid point to the difference equation 
at that same point. That is to say, we change that unknown 
to satisfy the corresponding equation; or, as in line 
relaxation, we change simultaneously a set of unknowns to 
satisfy the corresponding set of equations. This ’Viarriage" 


between the unknown and the equation at the same point is 
not always natural. In many cases, especially in solving 
a system of differential equations (q >1), the natural 
thing is to change several unknowns in ord; r to satisfy just 
one difference equation. Such a scheme is called distributed 
relaxation. (See Lecture 7 in Brandt (1978a).) A special 
case of such a relaxation was suggested by Kaezmarz (1937) 

9) 

and analyzed by Tanabe (1971) . Various cases of distrib- 

uted relaxation for singular-perturbation problems are 
analyzed by Dinar (1978). 

Let us show how distributed-relaxation yields uniformly 
bounded smoothing rates even when the marching direction 
is upstream, i.e., against the direction of convection. 

Take again the operator (6.3a) and assume forward relaxation. 
Instead of changing only the approximation u^(x) tc 
satisfy L^u^‘(x) = T(x), change now both u^(x) and 
u^(x + h) : Change u^(x) by adding to it 6, and 
u^(x + h) by adding to it -w5, where w is a fixed 
coefficient (see below) and 5 is calculated so that the 
equation L^u^(x) = F(x) is satisfied after these changes. - 
This marching process is stable for w < (1 + q)/2. The 
larger w the better is the smoothing rate. By taking w 
not far from the critical value (1 + n)/2, we can get 
smoothing rates v which are less than 1 for all q, and 
V = 0(ri ^) for large n. 

w is called the distribution coefficient, and should not 
be confused with the familiar "relaxation parameter". The 
above scheme is called Distributed Gauss-Seidel (DCS) 
relaxation, because, as in the Gauss-Seidel scheme, each 
difference equation in its turn is fully satisfied by the 


changes. For all problems examined, 'ncluding incompressible 
Kavier-Stokes equations, DCS was found to be the best 
smoother. 


One final remark: The difference operat.or must be 

uniformly stable (see Section 5) , otherwise no relaxation 
scheme can have uniformly bounded smoothing rates. For 
example, if the central difference approximation is used 
instead of (6.3a), ev'en backward relaxation would have 
1 - n/2 


u (n) = 


2 + i(l + n/2) 


-^1 as 


FOOTNOTES 

exception is the case when the coarse-grid difference 
H 

operator L does not fully use the smoothness of the solu- 
tion. In that case, if I, in (2.4) is of sufficiently 
high order, then V will be smooth enough to be appfoxi- 
mated by some V . This situation is related, however, to 
the use of an inappropriate approximation order, and w’ill 
therefore not arise in a fully adaptive procedure. 

Provided the tv?o appearing in (2.9) are 

identically the same . A coinaon programming error is that 
they differ at some special points. 

3) 

It is not necessary to compute the residual norm, since 
this particular algorithm is "fixed", its flow does not 
depend on the internal measures, and the number of sw’eeps 
made at each stage is prescribed in advance. For more 
complicated equations an "accommodative" algorithm, with 
internal switching criteria (e.g., the algorithm in Figure 1 
above), may be desired. But, for more complicated equations, 
relaxation is more expensive, so that the extra w’ork in 
computing [[r^jj is relatively small. 


^^x-cxtrapolalion is produced by the saaie FASF^IG program 
of HUGTArE (1978) through simple changes shovni there by 
Comment cards. 

^^Alternati . ely , extra smoothness on the scale of the 
finest grid can be used to produce a solution with 

errors smaller than the truncation errors in very little 
work. Indeed, if the difference equations do not exploit 
all the smoothness in the solution, an approximation to the 
level of the G^ truncation errors is obtained (with 
T-extrapolation) already on one of the coarser grids. All 
is needed then is to interpolate from that grid to G , 
with high enough order of interpolation. 

^^For example, in the approach taken by Hackbusch (1978), 
the solution of a coupled pair of elliptic equations ' . 
requires v’ork equivalent to too many (at least 28/3, 
instead of just 2) solutions of a single equation. 

^^An abnormal run can usually be detected by examining 
the condensed output (output similar to the first three 
columns in Table 1). See Debugging Techniques, Lecture 18 
in Brandt (1978a). 

is enough to study relaxation schemes for the lowest 
order operator, because (i) one can compute higher— order 
approximations via the lower-order ones (see Section 3.4). 
(ii) For any relaxation scheme, the smoothing-rate depen- 
dence on the approximation-order is not very significant. 

References due to Blair Swartz and Gene Golub. 

10)ph definition of r'' should later be under- 
stood as the current right-hand side f^ «= ^h/2 Fig. 

1 ).‘ 
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